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1.  Introduction 

One  of  the  main  tasks  of  functional  genomics  is  the  study  of  the  role/activity  of  a  given 
expressed/purified  protein  in  vivo,  for  which  it  should  be  inserted  in  cultured  cells  or  an 
organism.  Initially,  this  problem  was  solved  in  an  indirect  way  based  on  the  transfection  of  cells 
with  a  proper  active  viral  DNA  (construction  of  a  proper  DNA-vector)  with  subsequent 
production  of  the  given  protein  by  cell's  own.  Later,  they  started  to  manage  without  the  DNA- 
transfection  stage,  trying  to  deliver  the  protein  into  the  cell  directly.  The  first  protein  delivery 
methods,  however,  were  the  same  as  those  for  DNA,  such  as  electroporation,  micro  injection, 
cationic  lipids,  and  others.  These  are,  in  particular,  too  invasive  and  low-efficient.  Thus, 
electroporation  often  entails  ruining  a  considerable  number  of  cells;  micro  injection  (insertion  of 
femtoliter  volumes  into  a  certain  cell  compartment)  needs  too  specific  equipment  and 
sophisticated  skill  while  only  single  cells  are  processed,  so  on.  In  addition,  these  methods  are 
suitable  for  laboratory  experiments  but  are  practically  inapplicable  in  vivo  for  their  potentially 
high  toxicity,  immunogenicity,  etc. 

That  is  why  in  the  last  decade  a  great  enthusiasm  was  sparked  to  the  discovered  ability  of 
some  proteins  (and  then  their  excised  sequences)  to  be  internalized  into  cells  by  themselves  and, 
moreover,  to  deliver  there  other  macromolecular  components,  often  much  larger  than  the  carrier. 
Early  reports  about  such  a  peptide  appeared  in  late  eighties  with  the  example  of  the  protein- 
transactivator  of  transcription  of  the  human  immunodeficiency  virus,  HIV-1  Tat  protein  [1-3].  In 
those  works  it  was  already  said  about  great  potentiality  of  the  phenomenon  for  the  intracellular 
delivery  of  macromolecules  previously  thought  to  be  impermeable  to  cellular  membranes;  the 
evidences  of  that  were  reported  soon,  see  e.g.  work  [4]  in  which  Fawell  et  al  showed  the  entry  of 
proteins  cross-linked  with  relatively  large  (at  that  time)  Tat-peptides  into  cultured  tissue  cells.  In 
a  sort,  a  boom  was  initiated  by  the  reports  of  Dowdy  et  al  in  late  nineties  (see  e.g.  [5])  about 
successful  and  fast  delivery  even  of  large  (e.g.  [3-galactozidase,  Mw  120  kDa)  proteins  fused 
with  a  small  peptide  (Tat-protein  fragment)  into  all  mouse  tissues,  brain  included,  after 
intraperitoneal  injection.  In  numerous  subsequent  works  one  can  trace  both  a  considerable 
extensive  progress  -  the  number  of  such  oligopeptides  called  protein  transduction  domains 
(PTDs)  has  grown  from  a  few  (HIV-1  Tat  PTD,  Drosophila  Antennapedia  homeodomain 
penetratin,  or  Antp,  and  some  others)  to  about  a  hundred,  as  is  noted  in  the  recent  mini-review 
[6]  -  and  a  diversity  of  serious  controversies,  presently  culminating  to  be  resolved.  The  reason  of 
such  a  sharp  interest  to  the  problem  is  obvious  at  least  from  its  applied  aspect.  From  a  number  of 
works  it  followed  that  the  translocation  of  such  peptides  and  their  cargo  seemed  practically  non- 


specific  and  almost  universal  (either  to  the  cell  type  or  cargo),  did  not  go  in  the  classical 
receptor/endocytosis  ways  and,  in  addition,  proceeded  much  faster  in  comparison  with 
conventionally  used  methods.  Also,  the  degree  of  risk  or  toxicity  turned  out  to  be  incomparably 
lower.  In  the  popular  scientific  literature  and  on  the  sites  of  commercial  producers  they  started  to 
speak  about  a  new  era  in  pharmacology  and  therapy  ([7-8],  see  also  [9]).  On  the  other  hand,  the 
phenomenon  looked  absolutely  incomprehensible  at  least  from  the  commonly  accepted 
physicochemical  views  on  the  impermeability  of  lipid  membranes  for  such  highly 
charged/hydrophilic  molecules  as  cell-penetrating  peptides  (CPPs). 1  Besides,  a  series  of  results 
have  been  called  into  question  for  unreliability/limitations  of  the  observation  methods,  especially 
those  that  denied  the  endocytic  pathways  of  the  translocation;  that  is,  called  into  question  was  the 
very  phenomenon  of  protein  transduction. 

The  goal  of  the  present  deliverable  is  to  review  the  current  literature  in  the  field  and  select 
those  of  the  data  available  that  can  serve  as  a  starting  point  in  the  physical  analysis  and  computer 
modeling  of  the  transduction  mechanisms  as  the  latter  aspects,  despite  of  the  problem 
importance,  still  remain  practically  inelaborate. 

2.  Cell-penetrating  peptides:  general  features  (outline) 

Langel  et  al  [6,11]  relate  the  tenn  CPP  to  up  to  thirty  amino  acid  amphiphilic  peptides 
which  can  be  internalized  by  cells  by  energy-independent  mechanisms  (although  the  endocytic 
pathway  is  not  exluded).  According  to  their  classification,  CPPs  can  be  conventionally  divided 
into  three  classes:  protein  derived  CPPs,  model  peptides  and  designed  peptides.  The  first  are  in 
fact  PTDs  (with  Tat  and  Antp  being  most  studied),  the  second  represent  sequences  mimicking 
the  structures  of  known  CPPs  (typical  are  oligoarginines,  e.g.  (R)n=7+9  [12]);  lastly,  the  third 
comprise  chimeras  composed  of  hydrophobic  and  hydrophilic  parts  of  different  origin  (e.g. 
transportan  consisting  of  a  part  of  galanin  attached  to  mastoparan  via  lysine). 

The  primary  structure  of  known  CPPS  exhibits  no  special  common  features  apart  from 
the  fact  that  all  of  them  are  positively  charged  and  amphipathic  (except  polycationic 
homopolymers).  Their  net  positive  charge  is  ensured  mainly  by  arginines  and,  to  lower  extent, 
lysines;  thus,  (+8)-charged  basic  domain  Tat49_57  RKKRRQRRR  contains  six  arginine  and  two 
lysine  residues. 

There  are  controversial  data  on  the  role  of  the  secondary  structure.  Thus,  Derossi  et  al 
[13]  revealed  that  internalization  of  Antp  is  not  hampered  by  introducing  up  to  three  prolines  into 


1  Most  general  tenn  comprising  both  PTDs  and  model  or  chimeric  oligopeptides  having  the 
mentioned  properties  [10]. 


2 


the  sequence,  i.e.  by  disruption  of  the  a-helical  structure.  Together  with  demonstration  of  the 
translocation  of  reverse  helices  and  D-enantiomers,  this  allowed  the  authors  to  conclude  that  the 
internalization  mechanism  is  receptor-independent  (but  see  [14]).  For  model  CPPs  studied  in 
[12]  (in  particular,  (Rfi  and  (rfi)  the  translocation  efficiency  is  somewhat  different,  but  the  main 
emphasis  was  put  on  the  guanidine  specificity  of  the  arginine  structure  and  the  length  and 
conformational  flexibility  of  side  chains  due  to  alkyl  spacers.  The  data  of  Zaro  and  Shen  [15] 
also  show  the  importance  of  guanidine  structure  for  transduction  while  mainly  the  number  of 
positive  charges  determines  endocytosis.  The  advantage  of  arginine-rich  peptides  over  lysine- 
rich  ones  and  minor  importance  of  other  structural  characteristics  are  rather  frequently  mentioned 
(see  e.g.  [16,17]),  though  there  are  some  reports  on  better  efficiency  of  lysine  homopolymers 
[18,19].  On  the  whole,  so  far  it  is  hard  to  come  to  a  more  or  less  definite  conclusion  on  the 
translocation  structural  preferences. 

Even  greater  diversity  of  data  can  be  reviewed  in  regard  to  the  types  of  cells,  unnatural 
lipid  vesicles  and  cargoes  transported  by  CPPs.  Overwhelming  majority  of  experiments  on 
successful  translocation  is  performed  in  vitro  on  cultured  mammalian  cells  grown  by  standard 
methods  and  then  incubated  with  solutions  of  CPPs.  As  for  cargoes,  their  spectrum  spreads  from 
small  peptides  and  oligonucleotides  (see  e.g.  [20])  to  proteins  of  120-150  kDa  [5,21]  and 
liposomes  of  hundreds  of  nanometers  [22]  (see  also  Tables  2,3  in  [6]).  They  can  be  attached  in 
different  ways  (usually  via  a  covalent  bond);  most  typical  are  shown  in  Fig.l  taken  from  paper 
[6]. 

Binding  and  translocation  mechanisms.  As  noted  above,  this  central  question  still 
remains  unclear  and  gives  rise  to  sharp  contradictions.  Starting  from  early  works  [2,24,25],  an 
intrigue  has  been  generated  by  reports  about  bypassing  classical  endocytic  pathways  and 
supposedly  direct  interaction  of  CPPs  with  the  lipid  bilayer  of  cellular  membranes.  The  main 
argument  in  favor  of  such  statement  was  the  preserved  ability  (of  either  TAT  PTD  or  Antp)  of 
translocation  at  low  (4  C)  temperature  and  also  under  removing  ATP  sources  of  energy  or  in  the 
presence  of  different  inhibitors  of  endocytosis.  Reports  of  this  kind  accumulated  till  2003 
[13,16,17,18,22,26],  involving,  apart  from  TAT  PTD  and  Antp,  the  growing  number  of  different 
CPPs,  although  opinions  in  favor  of  endocytosis  were  also  expressed  (e.g.  [27]). 

In  2003,  however,  the  situation  seemed  to  be  changed  due  to  appearance  of  several  works 
(first  of  all  [28])  pointed  at  serious  imperfections  of  the  most  widespread  methods  (cell  fixation, 
flow  cytometry,  etc)  resulting  in  the  fact  that  CPPs,  strongly  interacting  with  the  cell  surface, 
were  not  removed  even  by  repeated  washings.  Consequently,  flow  cytometry  does  not 


2  However,  the  number  of  arginines,  as  noted  in  [16],  should  not  be  excessive;  thus,  for 
polyarginines  n  =  8  is  opimal,  with  internalization  being  practically  abolished  at  n  >  12  . 
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discriminate  between  internalized  and  membrane-associated  peptides,  and  this  of  course  leads  to 
overrated  indices  of  penetration  and  other  artefacts  and  wrong  conclusions.  The  authors  of  paper 
[28]  insisted  that  in  live,  non-fixed  cells  a  typical  endosomal  distribution  of  peptides  could  only 
be  observed  and  the  uptake  is  inhibited  at  low  temperature  and/or  ATP  depletion.  Similar  results 
have  been  obtained  in  a  series  of  other  works  of  that  year  [29-32]. 

H 


Fig.l.  Attachment  of  cargoes  (green)  to  CPPs  (blue).  A,B,C  and  D  represent  cargo  covalently 
bound  to  the  CPP  via  a  peptide  bond,  thiazolidine  ring,  disulphide  bridge  and  bifunctional  linker 
molecule,  respectively.  ( E)  A  large  cargo  molecule  (e.g.  streptavidin,  shown  in  dark  green)  is 
non-co valently  bound  to  a  smaller  cargo  (e.g.  biotin,  light  green)  that  is  covalently  attached  to 
the  CPP  [6], 

Meanwhile,  in  the  same  year  at  least  two  papers  appeared  [15,39]  whose  authors  took  the 
criticism  expressed  in  [28-32]  into  account  and,  nevertheless,  came  to  an  unambiguous 
conclusion  in  favor  of  the  existence  of  an  energy-independent  non-endocytic  channel  of 
internalization  at  least  for  several  CPPs  (e.g.  (R)7),  though  for  TAT  PTD  and  Antp  this  channel 
was  found  to  be  somewhat  restricted.  In  work  [15]  with  the  help  of  a  newly  proposed  method  of 
subcellular  fractionation  the  contributions  of  transduction  and  endocytosis  were  separated,  with 
the  transduction  ability  being  observed  for  different  polyarginine  peptides  and  Tat47.57.  As 
improbable,  the  endocytic  pathway  is  also  noted  in  works  [34,35]. 

To  the  moment,  one  can  register  a  temporal  dynamical  equilibrium  of  opinions  on  this 
central  question.  As  follows  from  recent  mini-review  [6],  the  transduction  ability  of  CPPs, 
although  remaining  incomprehensible  so  far,  is  doubtless  -  the  authors  of  [6]  in  fact  include  it 
into  the  very  definition  of  CPPs.  On  the  other  hand,  for  many  CPPs  the  role  of  endocytosis  is  far 
from  negligible  -  this  is  being  confirmed  by  new  evidences  [36-39],  -  though  the  endocytic 
modes  could  be  non-classical  (e.g.  caveolar,  lipid  raft,  macropinocytosis,  so  on  [29,40]). 
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There  are  relatively  less  data  on  more  detailed  models  of  CPP  membrane  binding  and 
translocation.  More  or  less  convincingly,  one  can  assert  that  an  important  role  in  the  CPP  binding 
to  the  cellular  surface  is  played  by  ubiquitous  proteoglycans  (PG),  precisely,  their 
glycosaminglycan  chains  (GAG),  most  frequently  heparan  sulphate  (HS).  Evidences  to  this, 
qualitative  estimations  of  the  binding/dissociation  constants  included,  are  given  in  many  works 
[18,19,27,29,41,42,43],  but  see  [26]  (there  also  exist  data  on  binding  of  e.g.  penetratin  to 
membranes  of  liposomes  [44]).  Here  the  main  role  belongs  to  electrostatic  interactions  and 
neutralization  of  opposite  charges  in  the  formed  complexes  of  positively  charged  PTDs  and 
negatively  charged  HSs.  In  early  works,  in  which  a  direct  interaction  of  PTD  with 
phospholipides  of  the  plasma  membrane  was  supposed,  the  inverted  micelle  mechanism  was 

suggested  [13],  see  Fig. 2.  This  model  gained 
reminiscences  in  recent  work  [45]  where,  in  a  word,  a 
reversing  of  the  membrane  (displacement  of 
phosphatidylserine  to  outer  surface  of  the  cell 
membrane)  resulted  from  TAT  PTD  transduction  was 
registered.  Allusions  to  possible  formation  of  a 
membrane  pore  are  relatively  rare  since  the  latter  of  a 
needed  size  presumably  would  be  fatal  for  the  cell. 


Fig.2.  The  peptide,  represented  as  a  dimer,  recruits 
negatively  charged  phospholipids  (fdled  circles)  and 
induces  the  formation  of  an  inverted  micelle.  The 
hydrophilic  cavity  of  the  micelle  accommodates  the 
peptide  and,  possibly,  sequences  attached  to  it  that  can 
subsequently  be  released  in  the  cytoplasmic 
compartment  [13]. 

Kinetic  analysis  of  the  translocation  process  is  hindered  by  experiment  difficulties,  that  is 
why  data  on  kinetics  are  scanty  and  dispersive.  The  general  conclusion  can  be  reduced  to  a 
relatively  fast  internalization  of  peptides  detected  in  the  plasma  membrane  and  cytosol 
sometimes  even  in  several  minutes,  and  in  2CH-60  min  -  often  in  the  nuclear  membrane/nucleus. 
To  the  moment,  as  is  justly  mentioned  by  Zorko  and  Fangel  [6],  it  is  reasonable  to  consider  a 
phenomenological  kinetic  scheme  like  shown  in  Fig. 3  with  all  its  stages  being  sufficiently 
justified.  From  their  point  of  view,  however,  even  this  scheme  is  too  complicated  and  needs 
further  simplification.  Naturally,  such  a  kinetic  analysis  could  yield  indicative  results  only  for 
revealing  the  translocation  mechanism.  But,  astonishingly,  even  this  necessary  step  is  practically 
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absent  in  the  whole  body  of  literature  in  the  field,  where  attention  is  paid  mostly  to  qualitative  or 
physiological/medical  aspects  of  experimental  results.  Although  it  goes  without  saying  that 
targeted  application  of  the  phenomenon  to  drug  delivery  and  clinical  therapy  will  be  impossible 
without  full  comprehension  of  its  physicochemical  mechanism. 


Fig.3.  Simplified  kinetic  scheme  for  CPP  internalization.  Indices  in  and  out  represent  the 
portions  of  CPP  or  its  degradation  products  inside  and  outside  the  cell,  respectively;  M  denotes 
membrane  bound  CPP;  Free  means  internalized  but  non-bound  CPP  (e.g.  in  cytosol);  Bound 
means  the  fraction  of  CPP  that  is  interacting  with  inner  cell  structures  (intracellular  membranes, 
proteins,  etc);  Degradation  products  result  from  proteolytic  cleavage  of  the  CPP  in  the  cell;  EC 
denotes  endocytosis  [6]. 

3.  Planned  ways  and  methods  of  further  investigation 

Resuming  the  previous  section,  one  can  note  the  following.  First  works  done  on  cell- 
penetrating  peptides  (CPPs)  can  be  classified  as  qualitative  observation  of  the  peptide  uptake  by 
cells.  Leaving  aside  the  questions  of  adequate  experimental  protocols  it  is  possible  to  postulate 
that  the  effective  uptake  of  pure  CPPs  or  different  cargoes  linked  to  CPPs  is  now  established 
with  no  doubt.  However,  apparent  uptake  of  CPPs  does  not  mean  that  the  peptides  cross  the 
membrane  in  a  “mysterious”  way  -  the  uptake  can  be  caused  by  the  well  known  endocytic 
pathway.  There  is  a  significant  body  of  studies  that  confirm  either  endocytic  internalization  or 
direct  membrane  crossing  by  CPPs,  but  they  can  be  hardly  considered  as  absolutely  convincing. 
It  was  postulated  recently  that  several  most  widely  studied  CPPs  such  as  penetratin,  TAT  and 
synthetic  oligoarginines  are  internalized  concurrently  by  both  endocytosis  and  direct  membrane 
crossing.  It  was  established  that  the  binding  of  CPPs  to  the  membranes  is  non-specific  and  is 
necessary  for  both  direct  crossing  and  endocytosis. 
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The  endocytic  uptake  of  CPPs  is  of  little  interest  from  the  point  of  view  of  the  current 
project.  In  this  case  the  peptide  resides  at  the  same  extracellular  compartment  and  the  very 
complex  energy-dependent  process  of  vesicle  formation  facilitates  the  uptake.  In  contrast,  direct 
translocation  of  CPPs  through  the  membrane  is  of  great  interest.  During  this  process  the  soluble 
CPP  molecule  manages  to  cross  the  hydrophobic  core  of  the  membrane.  This  requires  some 
unknown  physical  mechanisms  making  it  possible  to  overcome  a  prohibitively  high  energy 
barrier  of  such  translocation.  That  is  why  we  henceforth  focus  on  direct  translocation  of  CPPs 
through  the  membrane. 

General  physical  picture  of  translocation.  In  general,  the  translocation  mechanism  can 
be  divided  into  several  stages  regardless  of  its  molecular  details.  The  first  stage  is  the  binding  of 
the  CPP  molecule  to  the  membrane  surface.  This  process  is  energetically  favorable.  The  next 
stage  is  translocation  through  the  hydrophobic  core  of  the  membrane.  This  stage  most  likely 
involves  the  crossing  of  a  substantial  energy  barrier  or  several  such  barriers.  After  this  stage,  the 
peptide  appears  bound  on  the  internal  side  of  the  membrane.  The  last  stage  is  the  release  of  the 
peptide  to  the  cytosol  that  is  again  energetically  unfavorable  because  of  the  binding  energy.  The 
translocation  is  spontaneous  and  does  not  require  external  energy  sources  like  ATP,  thus  the 
whole  process  can  be  viewed  as  a  thermo-activated  diffusion  process.  It  is  possible  to  introduce 
the  translocation  coordinate  x,  which  can  be  roughly  associated  with  the  position  of  the  CPP  in 
the  membrane.  A  hypothetic  energy  profile  along  the  translocation  coordinate  consists  of  two 
energy  minima  that  correspond  to  the  CPP  bound  to  the  external  or  internal  sides  of  the 
membrane.  These  minima  are  divided  by  the  energy  barrier  that  corresponds  to  the  peptide 
location  in  the  hydrophobic  core  of  the  membrane.  The  exact  shape  of  the  energy  profile  in  this 
region  depends  strongly  on  molecular  details  of  translocation.  Several  possible  variants  are 
discussed  below. 


Fig.4. 


As  noted  above,  there  are  several  models  of  translocation  mechanisms  in  the  literature 
(inverted  micelle  model,  carpet  model,  etc).  It  is  also  possible  that  the  single  peptide  intercalates 
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between  the  lipids  and  cross  the  membrane  without  forming  specific  pore-like  structures.  Let  us 
consider  these  cases  in  more  details. 

Inverted  micelle  model  and  carpet  model.  These  imply  that  the  lipid  bilayer  is 
reorganized  in  such  a  way  that  no  contact  between  the  peptide  and  the  hydrophobic  lipid  tails  is 
necessary  during  translocation.  In  fact  the  peptide  is  coated  by  the  lipid  heads. 

The  inverted  micelle  model  implies  that  the  peptide  molecule  induces  the  formation  of  a 
local  depression  on  the  surface  of  the  membrane.  This  can  be  accomplished  by  strong  binding  of 
the  peptide  to  the  head  groups  of  several  lipids.  The  depression  can  then  close  itself  into  an 
inverted  micelle.  If  this  micelle  is  opened  to  the  internal  solution,  the  peptide  is  transferred 
through  the  membrane.  However,  the  feasibility  of  this  mechanism  is  questionable:  the  minimal 
size  of  inverted  micelle  is  quite  large,  leading  to  severe  distortion  of  the  outer  leaflets  of  the 
bilayer;  thus  the  structure  can  hardly  be  stable. 

The  carpet  model  implies  that  a  number  of  peptide  molecules  cover  the  surface  of  the 
membrane  as  a  patch  of  a  carpet.  Interactions  between  the  peptides  and  lipids  destabilize  the 
planar  structure  of  the  bilayer  and  induce  the  fonnation  of  a  pore,  covered  by  the  lipid  molecules 
and  peptides.  Such  a  pore  does  not  expose  hydrophobic  lipid  tails  to  the  water  environment  -  it  is 
covered  by  the  lipid  head  groups  with  bound  peptides. 

The  barrel-stave  model,  in  contrast,  postulates  that  the  peptides  themselves  form  the  pore 
walls.  The  leaflets  of  the  lipid  bilayer  in  this  case  are  not  bent  at  all.  The  peptide  molecules  are 
arranged  like  a  barrel  with  the  hydrophobic  outer  surface  contacting  the  lipids  and  the 
hydrophilic  inner  surface  that  forms  the  pore  [46]. 

All  the  mentioned  structures  should  be  meta-stable  with  the  lifetime  comparable  to  the 
translocation  time.  However,  fonnation  of  such  structures  is  energetically  unfavorable  because  it 
requires  a  significant  perturbation  of  the  bilayer.  As  a  result,  the  energy  profile  for  these  models 
should  be  like  that  in  Fig. 5: 
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Direct  translocation.  If  the  peptide  intercalates  between  the  hydrophobic  tails  of  the 
lipids  and  crosses  the  membrane  directly  without  formation  of  a  specific  membrane  structure 
then  the  energy  of  peptide  inside  the  membrane  core  should  be  the  highest  and  the  energy  profile 
will  have  a  single  barrier  (Fig. 6). 

Peptide  between 


Quantitative  parameters  of  binding.  Binding  of  the  CPP  to  the  membrane  surface  is  the 
crucial  stage  of  translocation,  presenting  in  any  translocation  model.  Recent  advances  in 
experimental  techniques  allow  one  to  detennine  the  binding  properties  quantitatively.  To  our 
knowledge,  most  of  recent  quantitative  data  relate  to  penetratin  or  Tat  PTD. 

The  binding  parameters  can  be  evaluated  using  the  binding  isotenns.  The  only  parameters 
which  can  be  directly  extracted  from  the  experiment  is  the  membrane-bound  peptide  to  lipid 
molar  ratio  and  apparent  dissociation  constant  of  the  peptide.  Using  the  model  of  Langmuir 
adsorption  it  is  possible  to  extract  the  binding  energy  E/,  by  fitting  the  binding  isotenns.  A  very 
interesting  attempt  to  extract  different  components  of  the  binding  energy  was  published  recently 
[44].  The  authors  use  the  Gouy-Chapman  theory  that  allows  them  to  estimate  the  electrostatic 
and  hydrophobic  components  of  the  binding  energy.  In  addition  it  is  possible  to  extract  such 
quantities  as  the  local  concentration  of  the  peptide  near  the  membrane,  membrane  surface  charge 
density,  membrane  surface  potential  and  effective  peptide  charge.  These  can  be  used  in 
modeling.  They  are  especially  useful  because  the  membrane  and  peptide  charges  can  be  adapted 
to  represent  not  only  the  outer  solution,  but  also  the  cell  interior,  which  is  necessary  to  calculate 
the  binding  to  the  cytosolic  side  of  the  membrane. 

Kinetics  of  internalization.  As  noted  above,  internalization  of  penetratin  and  other  CPPs 
includes  two  mechanisms  -  direct  penetration  and  endocytosis.  Zaro  and  Shen  [15]  published  the 
experimental  data  that  allows  one  to  compare  the  importance  of  these  pathways  quantitatively. 
They  measured  the  time-resolved  concentration  of  the  internalized  peptide  in  the  vesicular 
fraction  and  in  the  cytosol  separately.  The  authors  did  not  calculate  any  kinetic  constants,  but 
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their  data  are  perfectly  suited  for  such  an  analysis.  The  only  serious  limitation  is  a  very  small 
number  of  points  on  the  curves,  which  can  make  the  fitting  ambiguous.  Let  us  consider  the  time- 
dependent  concentrations  of  the  peptide  in  the  vesicles,  V(t),  and  in  the  cytosol,  C(7),  measurable 
experimentally.  The  peptide  concentration  in  the  outer  solution  is  O  (Fig.  7). 


Fig.7. 

The  processes  of  direct  penetration  and  vesicle  formation  can  be  approximated  by  single- 
stage  first-order  reactions.  One  can  define  the  kinetic  constants  of  peptide  exchange  between  the 
three  compartments  (cytosol,  vesicles  and  outer  solution),  kcv,  kvc,  kco,  koc,  kvo,  kov. 

The  following  balance  equations  describe  the  behavior  of  the  system: 

C  =  kvcV  +  kocO-(kcv  +  kco)C 
V  =  kcvC  +  kmO  -  (kvc  +  kvo)V 

This  set  contains  six  empirical  kinetic  constants.  The  quality  of  experimental  data  does 
not  allow  us  to  determine  all  the  six  constants  by  direct  fitting  the  curves.  That  is  why  further 
simplifications  are  needed.  The  endocytosis  is  an  energy-dependent  process  that  can  be 
considered  as  irreversible,  thus  the  vesicles  can  hardly  bring  their  cargo  back  to  the  outer 
solutions  and  we  can  set  kvo  =  0.  Next,  it  is  possible  to  assume  that  the  direct  permeation  through 
the  membrane  of  the  vesicle  is  essentially  the  same  as  through  the  outer  membrane  (in  reality  this 
can  be  modified  by  the  acidic  conditions  inside  the  vesicles).  In  this  case  kvc  =  koc  =  kin;  kcv  =  kco 
=  kout.  The  remaining  three  constants  can  be  detennined  by  fitting  the  experimental  data. 

Determining  the  energy  profile  of  direct  translocation.  Analysis  of  the  kinetics  allows 
one  to  obtain  reliable  quantitative  kinetic  constants  k„  and  kout  that  describe  the  inward  and 
outward  penetration  through  the  membrane.  According  to  the  Kramers  theory  of  thermo- 
activated  reactions,  one  can  write: 
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K  =  A  exp(-  Ein  /  kgT) 

kout=  A^v(~E out1  kBT) 


where  Ein  and  Eout  are  the  heights  of  the  effective  energy  barriers  of  penetration,  A  is  the  Kramers 
factor.  The  effective  energy  barrier  depends  on  the  heights  of  actual  barriers  along  the 
translocation  coordinate.  If  there  are  no  significantly  populated  intermediate  states  on  the 
translocation  pathway,  then  the  effective  barrier  corresponds  to  the  rate-detennining  highest 
barrier  along  the  translocation  coordinate. 

The  heights  of  the  other  energy  barriers  cannot  be  extracted  using  kinetic  data  only.  In 
order  to  do  this  the  whole  body  of  experimental  data  should  be  involved.  The  following  scheme 
(Fig.  8)  shows  which  kind  of  data  can  be  used  to  estimate  the  heights  of  various  barriers: 


Fig.8. 


Theoretical  estimate  of  the  rate-determining  energy  barrier.  The  values  of  Ein  and  Eout 
that  are  most  probably  the  rate-detennining  barriers  heights  obtained  from  experimental  data 
become  invalid  if  different  peptides  or  different  membrane  composition  are  considered.  Thus, 
independent  theoretical  estimates  are  necessary.  One  of  the  possible  approaches  is  direct 
molecular  dynamics  calculations  of  the  molecular  membrane  models  corresponding  to  various 
penetration  mechanisms.  This  approach  can  discriminate  between  possible  mechanisms  and 
determine  the  most  probable  mechanism  for  a  particular  membrane  composition,  penetrating 
peptide  and  cargo.  However,  such  calculations  require  the  modeling  of  quite  large  membrane 
patches  for  quite  long  periods  of  time  that  presumes  the  usage  of  supercomputers  containing  not 
less  than  100  processors  for  several  months. 

Alternatively,  simplified  theoretical  models  can  be  developed  for  each  penetration 
mechanism.  This  approach  does  not  require  extensive  calculations  but  is  error-prone  because  of  a 
number  of  empirical  data  incorporated  into  such  models.  The  lack  of  quantitative  data  does  not 
allow  one  to  extract  the  empirical  constants  with  proper  accuracy  for  a  wide  range  of  membranes 
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and  penetrating  peptides.  Further  progress  needs  experiments  closely  coupled  with  theoretical 
studies  and  directly  targeted  on  the  verification  of  specific  theoretical  models. 

General  approach  and  formalism.  To  construct  a  quantitative  model  of  the  phenomenon 
we  will  consider  the  translocation  process  as  a  diffusion  process  with  allowance  for  structural 
properties  of  the  (CPP  +  cargo)  complex,  its  interaction  with  the  membrane,  and  the  structural 
and  dynamical  properties  of  the  membrane.  Here  we  outline  the  possible  way  of  perfonning  this 
task. 

The  diffusion  process  of  crossing  the  membrane  by  the  CPP  is  conditioned  by 
fluctuations  in  positions  of  structural  components  of  the  CPP  and  membrane.  This  implies  the 
two  most  applicable  methods  of  description:  Langevin  stochastic  equations  or  corresponding 
Fokker-Planck  equations  for  the  multidimensional  distribution  function  of  system’s  variables.  Let 
Fcpp(x)  be  potential  expressed  in  tenns  of  the  CPP's  structural  variables  x  =  (x1,x2...xAf)  which 

describe  the  relative  motion  of  CPP's  structural  components  and  motion  of  the  CPP  as  a  whole. 
Next,  let  VM(x)  be  potential  detennining  the  dynamics  of  structural  components 

y  =  y2...  >'w)  of  the  membrane,  and  W(x,y)  stands  for  the  interaction  between  the  CPP  and 

membrane.  Then  the  Langevin  stochastic  equations  describing  both  dynamical  and  diffusive 
processes  in  the  system  read: 

^  +  pkBTq,m.  (mtjC'j)  =  S„S(t  -  f) 

=  -  +  V2  i=  {1...JV},  j  =  { 


where  U(x,y)  =  VCPP(x)  +  VM(y)  +  W(x,y)  is  the  total  potential  energy  of  the  system  (CPP  + 
M)  and  qf  is  the  friction  coefficient  for  z-th  coordinate. 

The  corresponding  Fokker-Planck  equation  for  the  joint  distribution  function  P(x,  y,t) 
reads  [47-48]: 


fs  a  r  1  8U(x,y)  |  d 

1 i=\  x‘  dxt  kBT  dx.  Ox, 


M 


a 

i 
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+a]1 
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where  the  diffusion  coefficients  Z> 

xi 


knT  J 
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relate  to  structural  variables  of  the 


CPP  and  membrane,  respectively. 
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In  both  cases,  in  view  of  the  real  number  of  variables,  to  directly  follow  the  equations 
above,  in  either  numerical  or  analytical  calculations,  is  a  daunting  task.  Here  extremely  helpful  in 
necessary  simplifications  are  the  methods  of  adiabatic  elimination  of  fast  variables  which  we 
successfully  applied  earlier  for  description  of  biomolecular  systems  [49-50].  We  also  intend 
using  the  Brownian  dynamics  methods  for  the  most  important  slow  system's  variables  (to  be 
developed  at  the  next  stage  of  the  project). 

Construction  of  a  dynamical  model  of  the  CPP  is  one  of  the  main  auxiliary  tasks.  Such  a 
model,  on  the  one  hand,  should  be  simple  enough  and  accessible  computationally.  On  the  other 
hand,  it  should  reflect  specific  features  of  a  given  CPP  that  determine  its  penetration  capability. 
Preliminary  analysis  shows  that  such  a  model  can  be  based  on  the  notion  of  the  CPP  as  a  system 
of  several  blocks  mutually  interacting  by  potential  Fcpp(jc).  The  relative  movements  of  the 
blocks  will  be  identified  by  degrees  of  freedom  x ,  and  the  corresponding  parameters  of  the  latter 
will  determine  the  object  specificity.  It  is  supposed  that  a  block  will  be  relatively  large  (at  a  level 
of  a  residue  or  larger)  to  decrease  the  number  of  degrees  of  freedom  to  reasonable  extent  but 
without  loss  of  the  CPP  specificity.  We  consider  the  use  of  the  coarse-grained  knowledge  based 
potentials  (see  e.g.  [511])  for  constructing  FCPP(jt)  as  acceptable  in  modeling. 

Not  less  important  is  the  way  of  modeling  the  membrane.  Here  we  also  intend  holding  the 
method  of  reasonable  sufficiency.  If  the  coarse-grained  model  is  exploited  in  modeling  the  CPP, 
then  the  continuum  approximation  will  be  sufficient  for  modeling  the  membrane  since  tracing 
the  behaviour  of  individual  lipid  molecules  is  not  necessary.  The  following  model  of  the 
membrane  can  be  suggested:  the  hydrophobic  central  part  of  the  membrane  (tails  of  lipid 
molecules)  is  represented  by  an  infinite  homogeneous  and  isotropic  dielectric  layer  of  finite 
thickness  dcore  with  a  small  dielectric  constant  ( s  =  2  -r-  3  )  taken  from  experiment.  This  layer  is 

covered  on  both  sides  with  two  dielectric  layers  of  thickness  d polar  that  correspond  to  the 

domains  of  polar  heads  of  lipid  molecules.  These  latter  layers  have  its  dielectric  constant  as 
typical  for  partially  ordered  polar  medium  (50  <£'<80).  The  polar  layers  possess  bulk  and/or 
surface  charges  representing  the  charges  of  real  lipid  heads,  The  charge  distribution  can  be  rather 
inhomogeneous  in  the  direction  normal  to  the  bilayer  but  isotropic  in  its  plane.  Such  a  model 
allows  us  to  calculate  the  potential  profiles  for  each  type  of  penetrating  particles  (structural 
blocks  inside  the  peptides).  It  is  possible  to  vary  the  lipid  composition  of  the  membrane  and 
asymmetry  of  the  internal  and  external  monolayers  by  means  of  varying  the  dielectric  properties 
and  charge  distributions  of  the  polar  layers.  At  the  same  time  the  model  admits  sufficiently 
simple  ways  of  taking  into  account  the  interaction  of  different  CPPs  with  the  membrane. 

The  simplest  can  be  carried  out  within  the  constant  potential  approximation, 
i9(x)  =  W(x,y0 )  .  Then  the  positions  of  the  membrane  structural  groups  are  supposed  to  be  fixed 
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at  given  mean  values  y0 ,  and  &(x)  determines  the  interaction  of  each  of  movable  elements  of 

the  CPP  with  the  membrane.  In  this  case  the  problem  is  reduced  to  the  study  of  the  over-crawling 
of  a  non-rigid  object  of  the  given  structure  containing  components  of  different  mobility  through  a 
rigid  pore  (cf.  e.g.  [52-54]).  More  sophisticated  description  allows  for  self-consistent  dynamics 
of  the  (CPP  +  membrane)  system.  It  is  the  version  of  the  theory  that  could  turn  out  to  be  most 
suitable  for  describing  the  striking  features  of  protein  transduction. 
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Summary 


During  the  last  decade  the  protein  transduction  domains  or,  more  generally,  cell- 
penetrating  peptides  (CPPs)  are  being  of  special  interest  as  a  promising  non-invasive  means 
of  delivery  of  macromolecular  components  and  whole  macromolecules  into  cells  and  cellular 
nuclei.  Although  extensive  experimental  data  have  been  reported  in  the  literature,  they  are 
highly  diverse  and  often  contradictory,  starting  from  the  very  possibility  of  CPPs’ 
transduction  and  especially  the  translocation  mechanisms.  We  review  the  current  state  of  the 
art  in  the  field  and  outline  possible  physical  approaches  to  the  modeling  and  quantitative 
analysis  of  the  phenomenon  of  transduction. 
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EOARD  Partner  Project  211  -  Deliverable  2 


1.  Introduction 

As  was  noted  in  Deliverable  1,  from  the  available  amount  of  data  it  is  hard  to  distinguish 
more  or  less  special  common  features  in  the  primary/secondary  structure  of  the  cell-penetrating 
peptides  (CPPs)  except  the  fact  that  all  of  them  are  positively  charged.  This  suggests  that  the 
modeling  of  the  CPP  at  the  specific  structure  level  can  be  preceded  by  a  consideration  in  which 
the  CPP  is  represented  by  a  simple  dielectrical  stereoobject  (sphere/ellipsoid,  cylinder,  etc)  with 
some  distribution  of  its  positive  charge.  The  electrostatic  interaction  with  the  negatively  charged 
surface  of  the  lipid  bilayer  membrane  will  naturally  course  approaching  of  the  CPP  and 
membrane,  accompanied  with  considerable  changes  in  the  structure  and  geometry  of  the  latter, 
up  to  formation  of  some  complex  whose  character  will  allow  us  to  evaluate  the  transduction 
possibility. 

Obviously,  in  such  problem  formulation  the  leading  part  is  played  by  the  method  of 
modeling  of  the  membrane  structure.  This  can  be  done  within  the  mechanical  continuous  models 
by  the  methods  of  thin  film  elasticity  theory,  etc  (see  e.g.  [1]).  The  main  parameters  entering 
these  models  are  material  constants  like  coefficients  of  surface  tension,  bending  elasticity, 
spontaneous  curvature,  so  on.  Although  such  approaches  made  a  good  showing  in  some 
problems  of  the  rearrangements  of  lipid  aggregates  (e.g.  bilayer-micelle  transitions),  it  should  be 
noted  that  rigorously  they  are  applicable  for  describing  small  membrane  deformations  only.  In 
the  presence  of  domains  of  essentially  different  lipid  packing  they  can  produce,  at  their  best,  only 
qualitatively  correct  results. 

Another  possibility  is  to  use  a  semi-microscopic  approach  based  on  the  so-called 
"molecular  lipid  models"  originating  from  works  of  Israelachvili  et  al  on  lipid  self-assembly  (see 
e.g.  [2-4]  and  "the  opposing  forces  model"  (OFM)  proposed  therein);  in  recent  years  it  is 
intensively  developed,  in  particular,  by  Ben-Shaul  and  co-authors  [5-9].  In  this  approach  the 
lipids  are  treated  as  some  stereoobjects  (cones,  cylinders,  etc)  able  in  certain  limits  to  change 
their  shape,  depending  on  the  type  of  packing  in  a  lipid  aggregate,  without  changing  their  volume 
(fluid  model  of  constant  density).  Within  such  models  all  the  main  forces  determining  the 
aggregate  shape  can  be  represented  by  a  small  number  of  relatively  simple  contributions  into  the 
free  energy  per  lipid  molecule. 

Thus,  the  modeling  strategy  at  the  present  stage  looks  as  follows.  The  free  energy  of  the 
system  (peptide  +  membrane)  is  minimized  by  variational  methods.  With  this,  both  the  free 
energy  profile  and  system  geometry  (changing  as  the  peptide  approaches  to  the  lipid  bilayer,  up 
to  critical  geometrical  changes  indicating  the  possibility  of  transduction)  can  be  found.  Typical 
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realistic  values  of  parameters  of  the  bilayer  representing  an  artificial  lipid  membrane,  as  well  as 
those  of  the  peptide,  are  more  or  less  known.  The  influence  of  the  electrolyte  surrounding  the 
membrane  and  peptide  can  be  considered  in  a  simple  Debye  model  or  in  a  more  detailed  way 
involving  the  Poisson-Boltzmann  equation  [7,  10-12]. 

Below  we  briefly  recall  the  necessary  knowledge  of  the  lipid  membranes  and  the 
opposing  forces  model  and  then  expound  our  modernization  of  this  model  and  the  method  of 
calculation  of  the  free  energy  which  will  be  exploited  at  the  next  stage  of  the  project. 

2.  Lipid  aggregates  and  the  extended  opposing  forces  model 

When  trying  to  create  molecular  models  of  lipid  aggregates  self-assembly  it  has  turned 
out  that  rather  good  results  (although  initially  at  a  qualitative  level  of  prediction  of  possible 
structures  only)  could  be  obtained  within  a  relatively  rough  phenomenological  approach  using 
mainly  the  lipid  geometrical  properties  only.  For  example,  one  can  compose  the  dimensionless 
packing  parameter  p  =  vl  a0lc  whose  value  immediately  indicates  the  type  of  a  possible  lipid 
aggregate. 


Fig.l.  The  surface  and  bulk  intermolecular  interactions  that  define  the  dimensionless  packing 
parameter,  or  average  molecular  "shape  factor"  p  of  amphiphiles  [4]. 

Here  v  is  the  volume  per  lipid  molecule,  lc  is  the  maximum  possible  extension  of  the 
flexible  hydrocarbon  chains,  and  a0  is  the  "optimum"  (i.e.  corresponding  to  minimum  free 

energy)  headgroup  area  at  the  hydrocarbon- water  interface,  see  Fig.  1 .  There  exists  the  so-called 
"generic  sequence"  of  the  main  amphiphilic  structures  and  corresponding  types  of  the  elementary 
lipid  volumes  appearing  with  p  growing: 
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with  numerous  more  complex  intennediate  structures  [4]. 

Stability  of  this  or  that  structure  is  detennined  by  a  subtle  balance  of  the  forces  on  the 
bilayer  surface  and  in  the  hydrophobic  core.  Nevertheless,  to  a  rather  good  extent,  this  balance 
can  be  represented  by  the  following  additive  contributions  into  the  free  energy  per  lipid,  /  .  The 
simplest  is  the  classical  expression  [2,3] 


/  -  fs+fh’ 


(1) 


describing  the  interplay  of  the  two  opposing  tendencies:  to  minimize  the  hydrocarbon-water 
interface  surface  due  to  hydrophobic  interactions  ( fs )  and  to  drive  the  hydrophilic  headgroups  to 

move  away  from  each  other,  maximizing  their  contact  with  water  ( fh  ).  In  the  opposing  forces 
model  the  contribution  fs  is  expressed  through  the  surface  tension  coefficient  ;/  : 


fs  =  Ya  i 


(2) 


where  a  is  the  area  per  lipid  headgroup,  and  fh  reads 


/»=-.  (3) 

a 

being  some  kind  of  treating  the  headgroup  interaction  within  a  two-dimensional  van  der  Waals 
equation  of  state. 1  The  competition  of  these  two  forces  results  in  the  existence  of  the  so-called 


1  The  repulsion  expressed  by  Eq.  (3)  can  reflect  not  only  the  elecrostatic  interaction  of  charged 
headgroups  but  also  other  effects  like  those  of  excluded  volume,  etc. 
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"optimum"  headgroup  area  a0  that  can  be  easily  found  by  minimization  of  /  (1):  a0  =  1 7  , 

so  that  for  small  deviations  of  a  from  a0  Eq.(l)  can  be  reduced  to 


f(a)  =  2ya0  +  —  (a-a0f.  (4) 

a0 

The  simplest  OFM  (1-4)  is  applicable,  rigorously  speaking,  to  lipid  self-assembly  into 
uniformly  packed  aggregates.  In  more  complex  cases  it  should  be  taken  into  account  that,  first, 
the  mentioned  opposing  forces  act  at  somewhat  different  surfaces  (see  Fig.l).  Accordingly,  there 
is  a  difference  between  the  headgroup  area  at  the  hydrocarbon-water  interface  af  and  that  related 

to  the  maximum  headgroups  interaction  surface  ah ,  with  the  distance  lh  between  these  surfaces, 
so  that 

/  =  /«;  +  —  •  (5) 

ah 

The  distinction  between  at ,  ah  and  a0  contributes,  in  particular,  to  the  lipid  layer 
curvature  effects.  Besides,  there  exists  the  conformational  contribution  of  the  lipid  tails  into  the 
free  energy,  fc  .  Given  non-compressibility  (constant  volume  per  lipid  v )  and  uniform  density  of 

"lipid  liquid"  in  the  hydrophobic  core,  this  contribution  can  be  often  written  as  fc  =  r(h-/c)2 

where  b  is  the  average  length  of  the  lipid  chain  (detailed  calculations  of  the  chain  packing  in  the 
statistical  mean-field  theory  confirm  the  validity  of  such  approximation  [9,13]).  Therefore,  the 
extended  OFM  is  represented  by  the  expression 

f  =  fs+  fh+  fc  =  Yai  +  —  +  r  (b  -  lc)2 .  (6) 

“h 

As  shown  in  e.g.  works  [6,8],  it  is  quite  applicable  to  the  description  of  rather  complex 
non-uniform  lipid  aggregates  and  lipid-protein  interactions. 

Typical  values  of  the  constants  entering  Eq.(6)  are  rather  well-defined  in  the  literature 
(see  e.g.  [4,8,14]).  Thus,  for  a  saturated  hydrocarbon  chain  lc  < /max  =  (0.154  +  0.126n)  nm,  and 

_ o  o 

v  =  (27.4  +  26.9 n)  -10"  nm  ,  where  /max  is  the  fully  extended  molecular  length  of  the  lipid  chain 

containing  n  alkyl  groups.  Some  of  these  values  used  in  calculations  of  amphiphilic 
hydrocarbon  systems  in  water  are  listed  in  Table  1. 
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Table  1 


Surface  tension  of  the  hydrocarbon-water 

interface  y 

0.12 kBT /  A1 ,  or  50mJ/m2  at  T  —  300  K 

Optimal  headgroup  area  a0 

50-10~19m2 

Lipid  volume  v  (a  lipid  with  two  Cis  chains) 

10-27m3 

Critical  chain  length  lc  («  =  12-rl8) 

10-h20A 

Packing  parameter  p  =  u/ a0lc  (n  =  18) 

0.8 

Lipid  tail  compressibility  r  (n  =  12  -r  16  ) 

0.11^0.08kgr/A2 

3.  Modification  of  the  extended  OFM 

For  the  process  under  study  that  suggests  considerable  changes  in  membrane  geometry 
we  modify  the  approach  described  in  work  [8]  and  applied  there  to  the  description  of  the 
structure  of  the  edge  of  a  planar  lipid  bilayer  and  calculation  of  the  corresponding  line  tension. 
Precisely,  in  the  approximation  used  below  we  assume  that  each  lipid  monolayer  of  the 
membrane  is  represented  by  a  layer  of  directed  molecules  of  two-dimensional  liquid  (smectic 
liquid  crystal).  The  length  and  tilt  of  the  molecule  director  determine  the  monolayer  thickness  at 
a  given  point.  So  far,  we  consider  only  one  type  of  lipids  composing  the  monolayers  and  neglect 
the  interaction  between  the  latter.  The  ends  of  the  lipid  tails  of  both  monolayers  are  situated  on 
the  membrane  "midline". 

Fig. 2  shows  a  schematic  model  of  such  a  bilayer  membrane.  According  to  supposed 
cylindrical  symmetry  of  the  problem,  we  imply  the  cylindrical  reference  frame  ( R ,  cp,  z)  with  the 
peptide  situated  on  the  symmetry  axis  z  .  In  this  frame  the  position  of  a  point  is  given  by  radius- 
vector  p  =  iR  cos  (p  +  jR  sin  <p  +  kz  where  i  ,j  ,k  are  the  orts  of  the  cartesian  frame  (x,  y,  z ) : 

x  =  Rcoscp 
y  =  f?sin<p 
z  =  z, 

where  R  =  x2  +  y1  is  the  distance  between  the  point  and  z  -axis,  and  <p  is  the  angle  between  the 
projection  of  the  radius-vector  p  onto  plane  (x,y)  and  x-axis. 
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We  now  need  a  proper  method  of  description  of  membrane  geometry  changing  as  the 
peptide  approaches  the  membrane.  As  distinct  from  that  used  in  work  [8],  we  take  into  account 
(i)  cylindrical  symmetry  of  the  system;  (ii)  nonequivalence  of  the  monolayers  and  complex  shape 
of  the  interlayer  surface;  (iii)  the  presence  of  the  peptide  of  a  given  geometry  and  charge 
distribution. 

3.1.  Geometrical  description.  As  a  vector  that  defines  the  average  position  of  the  inner 
end  of  a  lipid  molecule  we  take  radius-vector  p  of  a  point  on  the  membrane  midline  (in  fact,  the 
set  of  such  points  defines  the  mentioned  midline).  In  the  outer  ("upper")  monolayer,  the  average 
position  of  the  molecule  with  its  inner  end  at  point  p  is  denoted  as  b(p)  and  the  position  of  the 
outer  end  of  the  hydrocarbon  tail  (near  the  headgroup)  -  as  r(p)  =  p  +  b(p) .  Finally,  the 
position  of  the  headgroup  center  is  b(p)(  1  +  lh  / b) ,  where  Ih  is  the  distance  between  the  end  of 

vector  b  and  the  center  of  its  headgroup.  In  the  inner  ("lower")  monolayer  the  corresponding 
vectors  are  q(p) ,  t{q)  =  p  +  q(p) ,  and  q(p)(  1  +  /,,  /  q) ,  respectively  (see  Fig.2).  Of  course,  it  is 
supposed  that,  in  view  of  the  system  symmetry,  any  physical  value  is  ^-independent, 
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A(cp,R,z)  =  A{R,z) ,  so  that  the  problem  can  be  reduced  to  a  two-dimensional  one  in  plane 
{R,z). 

Suppose  that  the  interlayer  midline  is  described  by  function  z  =  z(R).  Let  us  parametrize 
all  the  necessary  functions  with  respect  to  R  as  a  parameter.  The  two-dimensional  radius-vector 
p  is  p  =  p(R)  =  eR  +  kz(R)  (where  e  =  R/R  =  i  coscp  +  jsirup  ),  and  d p(R)  =  (e  +  z' (R)k^ dR  . 
Here  the  prime  denotes  the  derivative  with  respect  to  R .  Consequently, 
b  =b(R)  =  ebR(R)  +  kb  (R)  where  bR  =b sin#  ,  b,  =b cos#  are  the  component  of  vector  b(R) 
and  6  is  the  angle  between  vector  b(R)  and  z  -axis  (tilt).  Then 


r(R)  =  p(R)  +  b(R)  =  (R  +  bR  )e  +  (z  +  b:  )k  , 


dr(R)=  (l-l -  b'R)e  +  (z' +  b'^k 


dR 


For  the  lipids  whose  inner  tails  fall  into  the  interval  ( R ,  R  +  dR),  the  infinitesimal  length 
of  the  line  of  the  hydrocarbon  interface  in  plane  ( R ,  z)  reads: 


dr\  =  dR^(l  +  bR)2  +(z'  +  b'z)2 


in  the  upper  layer  and 


di\  =  dR^{\  +  q'R)2  +(z'  +  q[f 


in  the  lower  layer.  Then  the  infinitesimal  area  on  the  upper  hydrocarbon  interface  is 


=  (R  +  bR)d(p\df  |  =  su{R)dcpdR  , 


where 


su(R)  =  (R  +  bR)J(l  +  b'R)2+(z'  +  b:)2  , 


and  on  the  lower  interface 


(V) 


(8) 
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dci'p  =  {R  +  qR)d(p\dt\  =  sl{R)d(pdR, 
where 

Sl(R)  ={R  +  ^r)^  +  ^'r)2  +(Z'  +  <l'z)2  ■ 


For  corresponding  infinitesimal  areas  on  the  upper  and  lower  headgroup  surfaces  we 
have,  respectively: 

da™  =  da™  ( R;b  +  lh ) ,  da f  =  dc a<°  (. R;q  +  lh ) .  (9) 


Proceed  to  the  infinitesimal  volume  of  the  membrane.  This  volume  occupies  a  cylindrical 
sector  resulted  from  rotation  of  the  above-considered  strip  in  plane  ( f?,z)  by  angle  d(p  around 
z  -axis.  Similarly  to  work  [8],  introduce  the  reference  frame 

Rs=R  +  &r{R) 
z^=z  +  %bz(R), 


where  ( kf ,  z^  )  are  the  coordinates  of  the  points  on  the  director  with  its  inner  end  localized  at 
( R ,  z) ;  variable  d,  takes  its  values  in  the  interval  (0,1)  (see  Fig. 3). 

For  the  upper  monolayer,  the  infinitesimal  volume  dV™  related  to  dg ,  as  is  clear  from 
Fig. 3,  can  be  written  as  dV™  =  R.d(p-bzdd,  ■  \dT% | ,  where  T\=p  +  %b,  so  that 


dr.  = 


(1  +  4b'R)e+{z'  +  b':)k 


dR  .  Therefore, 


dr.  =  b,dcpdR\ d4(R  +  4b, )  (b'R2  +  b'} )  +  2£ (b'K  +  z% )  + 1  +  z'1 

0 
1 

=  bJtpdR J  d%(R  +  %bK)  ^2A2+£Al+A0 
0 


where  A0  =  1  +  z'2 ,  =  2  [b'R  +  z'b[ ) ,  and  A,  =  b'R2  +  b'2 . 

Introduce  the  designations: 
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(24  4  )  yf^2  +  4  +  4  4  74 

44 


/(*)  =  j^7£24+£4  +  4  = 

o 


444-4 


x2^0 


84 


3/2 


In 


24  +  4  +  274(4  +  4  +4) 


4  +  2^44 


/1(;?)  =  }j4!7£24+£4+4 

0 


V(4  +  4  +  4)  --y/4 

34 


A 

24 


/(*)• 


Fig.3. 


Then 


JFU  =  bzdcpdR 


R-b„ 


A 


2A 


^  7(4  +4  +  4 )  —  7  4 

/(/?)  +  4  - i - ^2 - 2L2L 


2  7 


34 


vu(R)dcpdR , 


where 


v„(4)  =  4 


4 


24 


/(i?)  +  4 


>/(4  +  4+  4)  _ 74 


■2  7 


34 


(10) 
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is  the  density  of  the  volume  with  respect  to  parameter  R .  Similar  expressions  for  the  lower 
monolayer  can  be  written  by  the  replacements  b  — »  q,  u—>l . 


3.2.  Free  energy  of  the  bilayer.  According  to  Eq.(6),  for  the  upper  monolayer 


B 


f„(R)  =  rai(R)  +  —-+r[b(R)-iy, 

a„(R) 


where  / ( R )  is  the  free  energy  per  lipid  whose  inner  tail  falls  into  the  interval  (R,  R  +  dR) .  Then 


the  quantity 


f(R) 


V 


is  the  free  energy  density  (as  mentioned  above,  it  is  supposed  that 


o  =  const ).  Therefore,  the  free  energy  of  this  layer  reads: 


Fm1  =  =-T  d<3fu(R)v.(R)dR  =  — }  fJR)vu(KW . 

u J  u  i  i  v  i 


0  0 


To  specify  fu(R)  in  this  formula,  we  need  to  find  the  expressions  for  quantities  a{f)  and 
a{u) .  For  this,  we  write 


<p+S(p  R+AR 

a\u\R)=  J  {  su(R)d<pdR*su(R)A<pAR 

(p  R 
(p+lA(p  R+AR 

v=  J  j  vu(R)d<pdR*vu(R)A<pAR 

(p  R 


(ID 


that  is  approximately  valid  if  A R  «  Rmm  ,  where  Rmin  is  the  minimal  local  curvature  radius  of  the 
monolayer  in  the  course  of  membrane  defonnation.  Note  that  for  a  spherical  micelle  or  semi- 
toroidal  edge  this  condition  does  not  hold  rigorously.  In  such  cases  Rmm  ~  b ,  being  I  (H20  A  for 

real  lipids.  At  the  same  time  the  average  inter-lipid  distance  is  also  ~10  A.  Then  more  accurate 
calculations  of  the  area  per  lipid,  including  integration  along  the  membrane  midline,  should  be 
performed. 

If,  nevertheless,  the  approximation  (1 1)  is  valid,  then,  obviously, 


,(«) 


(R)  = 


Su(R>b)u 

\(R) 


and 


su(R-,b  +  Ih)^ 
\(R) 
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where  (R;  b  +  lh )  =  [R  +  bR  (l  +  lh  I  bj\  ^(l  +  b'R)2  +(z'  +  K)2  . 


Consequently,  the  free  energy  per  lipid  is 


fu(R)-rv^jj^+~-  +  T[b(R)-ic f 

vu(R)  V  su(R;b  +  lh) 


and  the  total  free  energy  of  the  upper  monolayer  takes  the  form 


Fm'  =~J/.(W  =^\\rs,(R,b)  +  ~  ,  s+lv'(R)[b(R)-lJ  \dR.  (12) 


v-  su(R;b  +  lh )  v 


For  a  planar  monolayer  the  free  energy  per  lipid  reads: 


B 


f0=ya-\ —  +  r(b-lcy  ,  ab  = 


v 


(13) 


(here  the  index  u  is  omitted  in  view  of  full  symmetry  of  a  planar  bilayer).  Also,  it  is  obvious  that 
dVu  =bRdRd(p,  da\u)  =  RdRdxp  ,  so  that  vu(R)  =  HR  and  su (R)  =  R .  In  this  case  /0  does  not 
depend  on  R  ,  and  the  free  energy  of  a  planar  monolayer  reads: 


77(°)  _ 

1  M 


(14) 


where  N  is  the  given  number  of  lipids.  The  value  b0  (or  a0  )  can  be  found  from  minimization  of 
/0  (13)  that  results  in  the  cubic  following  equation: 

ya3 -{B -2zuIc)a-2TV2  =0 .  '  (15) 

Obviously,  L  =  ^Nol  nbQ  =  yjNa0  /  n  . 

Introduce  dimensionless  quantities,  with  b0  serving  as  a  length  unit: 


t-L- 


b-B 


yu 


B- 


yv 
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Then  the  minimum  free  energy  per  lipid  in  the  planar  monolayer  (at  N  =  const )  is 


/•(min) 

J  0 


=  7 


o 

bn  L 


l  +  5  +  r(l-4)2 


Parameters  r ,  B,  lc ,  as  is  seen  from  Eq.(15)  at  a  =  a0,  are  not  independent  of  each  other, 
precisely, 

_  1  1  -B 


so  that 


/’(min) 
/  0 


V 


1 +  B  +  - 
2 


(16) 


Subtracting  now  (14)  with  /0  as  (16)  from  (12)  yields  the  free  energy  of  the  upper 
monolayer  deformation  A/7))' 1 .  It  could  be  also  written  in  dimensionless  quantities,  introducing 

R  =  R/ b0,  L  =  L/ b0,a  =  a  /  etc;  then  Ju(R)  =  su(R)/b0,  vu(R)  =  vu(R)  /  b^  where  JU(R)  are 
given  by  Eqs.(8),(10)  with  all  the  quantities  are  supplied  with  the  bar.  Finally, 


J  " 

o  l 

_su(R;b+lh) 

1  1  -B 


2  1  -/ 


(6TO-i)=v„TO-(i-4)A 


—  \2  - 


Naturally,  the  expression  for  the  free  energy  of  the  lower  monolayer  can  be  written  with 
replacing  u  /,  b  — >  q  . 


4.  Contribution  of  the  CPP  electrostatic  potential 

The  local  charge  density  on  the  maximum  headgroups  interaction  surface  of  the  upper 
monolayer  can  be  written  as 


12 


z 


°uW  = 


z 

ah(R) 


aj(R',b  +  lh) 


Zvu(R) 

su(R;b  +  lh)u’ 


where  Z  is  the  total  headgroup  charge,  and  the  other  quantities  are  found  above,  see  (7-10). 
Then  the  electrostatic  part  of  the  free  energy  originated  from  the  interaction  between  the  peptide 
and  the  upper  monolayer  reads: 


K 


ZZYdal 

ai(R',b  +  lh) 


2KZ'^Y(R)su(R:b) 
o  l  su(R;b  +  lh) 


vu(R)dR , 


(13) 


where  'F(r)  is  the  electric  potential  at  point  r  . 

The  real  distribution  of  the  electric  field  created  by  the  peptide  near  the  membrane 
depends  on  polarizability  and  composition  of  the  solvent.  The  usual  tools  in  allowing  for  these 
factors  are  mean-field  theories,  like  based  on  the  Poisson-Boltzmann  equation  or  its  simplest 
version  (Debye  screening).  At  the  initial  stage  of  calculations  we  restrict  ourselves  with  the 
latter. 

While  potential  'F(r)  ensures  attraction  between  the  peptide  and  charged  membrane 
surface,  the  peptide  surface  is  supposed  impermeably  rigid.  We  will  model  the  peptide  by  a  two- 
axis  ellipsoid  resulted  from  rotation  of  the  ellipse  with  axes  dxv  =  dR  and  <7,  around  z  -axis 

(with  the  ellipse  center  located  at  point  zp  and  the  ellipse  axis  <7_  running  along  z  -axis;  for 
definiteness,  dz  >  dR  ).  Then  the  ellipsoid  surface  is  composed  of  points  (x,y,z)  satisfying  the 
equation 

*2  +  v2  ,  (^-V>2_-, 
d\  d\ 


With  the  relationships  dz=sd,  dR  =  | yjs2  -1  j  d  ,  where  d  =  ^d2R  +dz  is  the  half¬ 
distance  between  the  ellipse  foci  and  s  =  ,  is  the  ellipse  elongation  parameter,  the 

4dl-d2R 

ellipse  surface  equation  can  be  re-written  as 

■t2  +  y2  ,  U-*P?  .2 

s1  e2-l 
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We  will  assume  that  there  are  n  elementary  charges  localized  at  points 

Qn  =  Xni  +  Ynj  +  Znk  {n  =  l,2,..,np^  inside  the  ellipsoid.  Then  the  electric  potential  at  point 

\ 

where  ad  is 

) 

the  corresponding  Debye  length. 

To  preserve  cylindrical  symmetry,  we  suppose  that  the  peptide  charges  are  uniformly 
situated  in  the  interval  {0,0,zp  +d  >  z  >  zp  -  dj  between  the  ellipse  foci  (Fig. 4).  Then 

r\  T 

Xn=Yn=  0 ,  and  Zn  =z  -d-\ - (n - 1)  .  Therefore, 

n  - 1 


exp 


where  rz(R)  =  z(R)  +  b„{R)  ■ 

For  the  membrane  surface,  obviously 

r(R)  =  /[/?  +  h/?(f?)]cos^9  +  j  [i?  +  hj?(f?)]sin^  +  k  [z(f?)  +  h-(f?)] 

=  e[R  +  bK(R)]  +  k[z(R)  +  b:(R)] 
r(R)  =  ^R2  +  b\R)  +  z2  (R)  +  bl+2  [f^(f?)  +  z(R)b:(R)]. 

The  proposed  way  of  defining  the  peptide  shape  and  charge  distribution  allows  us  to  vary 
the  shape  from  spherical  ( d  =  0)  to  cylindrical  (s  =  1)  as  well  as  the  distance  d(s- 1)  from  the 

last  charge  (located  in  the  focus)  to  the  lowest  ellipse  point  {0,0,z^ -dzj .  This  is  helpful  in 

modeling  real  peptides  of  different  size,  shape  and  charge. 


'F(r)  =  vF(r,rz)  =  J 


-  —  Jr2 -2rZ  +Z2n 

n  V  z  n  n 

_ 

Jr-2r:Z„+Z; 


=  rxi  +ryj  +  r.k,  created  by  these  charges,  is  '¥(r)  =  Vj — ^-^r-,exp 

n  V  R 


r-R 


Zn 
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z 


5.  Methods  of  computation 

We  now  need  to  obtain  the  system  free  energy  profile  as  the  peptide  approaches  the 
membrane.  This  task  will  be  done  by  minimization  of  the  free  energy  functional 
A Fu  +  Fmt  =  F(z(R),bR(R),bz(R),qR(R),q.(R))  with  respect  to  five  functions  z(R),  bR(R),  bz(R), 

qR(R),  and  qz(R)  for  every  fixed  position  of  the  peptide.  Obviously,  this  cannot  be  done 
analytically,  and  the  following  standard  numerical  will  be  exploited. 

The  independent  variable  R  is  discretized  as  Ru) ,  1  =  0,1. ..TV,  R{,)  =0,  R(N)  =  L.  The 
boundary  distance  L  is  taken  large  enough  to  ensure  that  all  functions  at  point  RN  are  the  same 
as  those  of  the  unperturbed  membrane.  This  results  in  the  following  set  of  boundary  conditions: 

z{L)  =  z\L)  =  0 ,  bR(L)  =  b'R(L )  =  b':(L)  =  0 ,  b£L)  =  b0 , 
qR(L)  =  q'R(L)  =  q[(L)  =  0  ,  qz(L)  =  -b0 , 

where  the  quantities  with  zero  subscript  are  related  to  the  unperturbed  membrane.  The  peptide 
center  is  located  at  point  ( 0,  zp ) .  All  the  functions  characterizing  the  membrane  are  properly 

discretized,  too.  The  i-th  discrete  part  of  the  membrane  is  characterized  by  five  parameters 
Zj,bRi,  ,  qRi,  and  qzi ,  with  each  being  varied  independently.  The  total  number  of  variable 
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parameters  is  5N  +  5 ,  and  functional  F  is  transformed  to  the  function  of  5N  +  5  variables, 
F(zw  ...z{N),bf0)  ...bfN),bz(0)  ...bz('N),qf0)  ...qR{N\qzw  ...qziN))  .  There  are  several  possible 

computational  approaches  to  this  multiparametric  optimization  problem.  We  will  consider  two  of 
them,  most  robust  and  easy  to  implement. 

1.  The  approach  of  local  variations  (also  known  as  unconditional  multi-parametric 
optimization).  Each  of  five  parameters  are  changed  by  small  quantity  ±A  at  particular  point  i 
and  the  new  value  Fnew  of  the  functional  is  computed.  If  Fnew  <  F ,  then  the  change  of  the 

parameter  is  accepted,  if  not  -  rejected.  The  procedure  is  applied  subsequently  to  all  N  points 
until  a  given  convergence  criterion  is  satisfied.  The  method  is  simple  in  implementation  but  very 
demanding  computationally  (although  several  techniques  are  known  to  improve  the 
performance).  The  boundary  conditions  are  easy  to  implement  in  this  method. 

2.  The  approach  of  non-local  variations  (shape  variations).  All  discretized  functions  are 
approximated  by  finite  sets  of  K  cosines 

z(,)  ='^jz(,)  cos(co(j)R{,}),  bfl)  =  's^bR,)  cos(coU)R(l)) ,  bU)  =  ^bzJ)  cos(co{J)R{,)) , 

7=0  7=0  7=0 

qR']  =  yifit(J)  cos  (cou>R{,)) ,  qf]  =  cos  (coU)R(,)) , 

7=0  7=0 

where  coi  is  a  set  of  "harmonics"  identical  for  all  parameters;  z( i] ,bf  '\q f'] , bz (j\qfn  are  the 
amplitudes,  of'1  =0.  As  a  result,  F  becomes  a  function  of  5tV  +  5  amplitudes 
F(z(0\..z(K),bR0\..bRK\bJ-°\..bz(K),qR0)  ...qRK),q_{0)  ...qz(K)) .  The  amplitudes  are  varied  to 
minimize  F  .  The  number  of  hannonics  can  be  quite  small.  The  reason  is  that  the  radius  of 
curvature  of  the  membrane  is  limited  -  it  is  larger  than  b0  but  smaller  than  L ,  thus  few 

harmonics  can  be  enough  to  approximate  any  possible  smooth  shape  of  the  membrane.  This 
means  that  the  number  of  variables  is  reduced  and  the  computational  intensity  is  substantially 
lower  than  for  the  scheme  of  local  variations.  In  addition,  all  derivatives  are  expressed 
analytically,  and  this  also  results  in  substantial  speed-up  of  calculations. 

The  boundary  conditions  are  transformed  to 

K 

£z0)  cos (tuu)L)  =  z(0) 

7=0 

K 

<  ’^jzU)cfj)  sin(c7(i)Z)  =  0 

7=0 
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The  remaining  eight  conditions  designated  by  dots  can  be  written  by  changing  z(j)  to 
bR(i)  ,qRU)  ,bzU)  and  qz]) . 

It  is  obvious  that  the  amplitudes  can  not  be  varied  independently  without  violating  the 
boundary  conditions.  At  least  two  amplitudes  should  be  changed  consistently  at  each  step.  This 
variational  scheme  is  much  more  complicated  in  implementation  as  compared  with  that  of  local 
variations.  The  usage  of  shape  variations  is  justified  if  the  computational  burden  of  local 
variations  scheme  becomes  intolerable. 
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Summary 


The  strategy  of  analytical  and  numerical  modeling  of  the  interaction  of  a  cell-penetrating 
peptide  with  a  lipid  bilayer  membrane  is  determined  and  substantiated.  As  a  basis,  we  choose  the 
extended  opposing  forces  model  which  we  essentially  modify  in  order  to  allow  for  considerable 
membrane  geometry  changes  induced  by  the  peptide.  The  membrane  is  considered  as  a  lipid 
fluid,  with  lipids  length  and  orientation  depending  on  their  positions.  Besides,  the  possibility  of 
deformation  of  the  membrane  as  a  whole  due  to  the  bending  of  the  interlayer  surface  is  also 
taken  into  consideration.  To  take  account  for  the  disturbing  influence  of  the  peptide  we  propose  a 
simple  three-parameter  model  of  the  latter,  allowing  us  to  vary  both  peptide  geometry  (that  can 
be  decisive  at  close  contact  with  the  membrane)  and  its  electrostatic  field. 

The  proposed  formalism  is  aimed  at  determination  of  the  system  free  energy  profile  and 
the  character  of  membrane  deformation  with  the  peptide  approaching.  This  will  be  done  by 
minimization  of  the  constructed  free  energy  functional  of  five  independent  functions  defining  the 
equilibrium  membrane  configuration  in  a  convenient  and  visual  way.  Several  minimization 
procedures  will  be  tested  to  find  the  optimal  one. 

We  suppose  that  the  free  energy  profile  to  be  found  will  allow  us  to  determine  the  optimal 
characteristics  of  the  peptide  in  regard  to  its  transduction  abilities.  This  will  also  make  it  possible 
to  construct  a  diffusion  kinetic  model  of  the  transduction  process  and  to  evaluate  its  kinetic 
parameters. 
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1.  Introduction 

Previously  (see  Deliverables  1,2)  we  have  analyzed  the  known  theoretical 
approaches  to  studying  the  dependence  of  the  shape  of  bilayer  phospholipid  membranes  on 
various  disturbing  factors  in  order  to  choose  the  most  promising  method  of  the  theoretical 
investigation  of  the  CPP  transduction  process.  As  a  result,  we  have  chosen  a  version  of  the 
phenomenological  opposing  forces  model  (OFM)  [1-3]  having  much  in  common  with  the  model 
of  oriented  liquid.  The  main  advantage  of  this  model  consists  in  its  comparative  simplicity  and 
supposed  ability  to  treat  considerable  deformations  of  the  bilayer.  True,  in  recent  papers  (see  e.g. 
[4])  the  OFM  was  used  mainly  for  description  of  small  perturbations  of  a  given  membrane 
geometry.  Needless  to  say,  studying  large  (even  if  local)  changes  of  the  membrane  shape,  one 
implies  that  the  model  exploited  is  able  to  describe  the  initial  unperturbed  membrane 
configuration.  At  the  same  time  such  a  model  must  be  instructive  enough  at  describing  the 
emergence  of  rather  complex  membrane  structures  of  nontrivial  geometry,  up  to  fonnation  of 
pores/channels  for  CPPs.  That  is  why  in  Delivery  2  we  have  presented  a  special  mathematical 
formalism,  aimed  directly  at  such  problems,  within  which  the  bilayer  lipid  membrane 
configurations  are  described  with  five  independent  functions  detennining  both  the  average 
characteristics  of  orientation  and  defonnation  of  the  lipids  in  each  monolayer  as  well  as  bends  of 
the  membrane  as  a  whole. 

The  present  stage  of  the  project  is  devoted  to  creating  the  programming  methods  of 
realization  of  the  model  developed  in  order  to  calculate  and  simulate  the  membrane 
configurations  at  various  values  of  the  model  parameters,  for  both  the  unperturbed  membrane 
and  that  interacting  with  the  CPP.  At  this  stage  the  problem  consists  in  numerical  minimization 
of  the  system  free  energy  functional,  with  the  membrane  configuration  given  by  the  five 
mentioned  functions  resulted  from  the  minimization  procedure. 

To  this  end,  we  go  two  different  ways.  The  first  is  to  apply  the  known  Euler-Lagrange 
formalism.  This  leads  to  the  necessity  of  numerically  solving  the  set  of  ten  nonlinear  first-order 
ODE  for  the  five  mentioned  functions  and  their  first  derivatives.  The  second  consists  in  direct 
variations  of  the  lipid  arrangement  in  the  membrane. 

2,  Minimization  of  the  membrane  free  energy  functional  by  the  Euler-Lagrange  method 

According  to  the  results  presented  in  Deliverable  2,  within  the  OFM  model  the  free 
energy  of  a  membrane  of  arbitrary  cylindric-symmetrical  configuration  reads 
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(1) 


F  = 
rM 


L 

(. R;z,z',bR,bz,b'R,b'z,qR,qz,q'R,q')dR 
o 


where  the  Lagrange  function  of  the  membrane 


L  (. R,z,z',bR,b'R,bZ’K,qR,q'R,qz,q'z)  =  Lu{R,z,z',bR,b'R,bz,b'z ) 

+  Lj  (R,z,z',qR,q'R,qz,q'z) 


(2) 


consists  of  the  two  tenns  related  to  the  "upper"  ( u )  and  "lower"  (/)  monolayer.  Functions 
bR(R),bz(R)  represent  the  corresponding  components  of  vector  b  ,  determining  the  length  and 
orientation  of  a  lipid  at  point  R  in  the  upper  monolayer  (see  Fig.2  of  Delivery  2),  and 
b'R(R),b'z(R)  are  their  derivatives  with  respect  to  R ;  qR(R),q,{R)  and  q'R(R),q'(R)  are  the 
analogues  for  the  lower  monolayer.  Lastly,  functions  z(R),z'(R )  detennine  the  shape  of  the 
midline  between  the  monolayers. 

Under  cylindric  symmetry,  for  e.g.  the  upper  monolayer  we  have  obtained: 


L  u(R,z,z',bR,b'R,b:,b':)  =  su(R;b)  +  B- 


v2u(R ) 


+  -•—  (b(.R)-lcfvu(R,b) 


su(R’b  +  lh)  2  l~l 

(see  Deliverable  2  for  the  conventional  constant  designations),  where 


(3) 


su(R,b)  =  (R  +  bR)yj(  1  +  b'R  Y  +  (z'  +  b’z)A  (4) 

with  su(R,b)dR  being  the  infinitesimal  hydrophobic  surface  corresponding  to  the  cylindric- 


symmetrical  infinitesimal  volume  vu(R)dR  ,  where 


v.(«) = b4lddR + 4bR)  V(i + f*;)2 + (*'+ (5) 

The  corresponding  analogues  for  the  lower  monolayer  can  be  written  by  the  replacement 
of  the  components  of  vector  b(R)  and  their  derivatives  by  those  of  vector  q(R) . 

The  Euler-Lagrange  equations  for  the  sought  functions  z(R),bR(R),b,(R),qR(R),qz(R) 

read: 
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These  quantities  expressed  in  terms  of  the  components  of  vectors  b(R)  or/and  q(R)  are 
given  in  Appendix.  Besides, 


=  — A)2  ~  y 4 ^23  —  y 6^25  ~  Tg^27  ~  Tl0^29 
B4=  L2-  L04  —  y4Lu  —  v6T45 
=  L06  -  y4LM)  -  v6T56 

B&=  L7  —  Z,og  —  v8T78  —  JioTgg 
-®10  =  Bg  -  L0l0  -  v8T710  -  v10T910. 


Inserting  (10)  into  (4)  yields 
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To  numerically  solve  set  (11)  we  have  developed  a  program  using  the  Runge-Kutta 
method  with  adaptive  step  size.  In  the  vicinity  of  the  zero  point  the  integrals  were  computed  by 
the  Gauss  method  with  prescribed  accuracy. 

We  have  obtained  preliminary  results  showing  that  there  exists  some  domain  in  the 
system  parameters  space  where  the  solution  can  be  found  in  a  narrow  interval  of  variable  R ; 
however,  extending  this  interval  leads  to  solution  instability  (see  Fig.  1  as  an  example).  The  lattr 
is  caused  by  a  series  of  singularities  like  zero  coefficients  at  the  highest  derivatives,  etc.  At 
present,  the  numerical  scheme  is  under  modernization  aimed  at  removing  such  drawbacks. 

3.  Minimization  of  the  free  energy  functional  by  the  method  of  local  variations 

At  present,  the  program  realizing  the  method  of  local  variations  (see  Delivery  2)  is  at  its 
testing-debugging  stage.  This  program  allows  us  to  obtain  the  optimal  configuration  of  the  free 
membrane  and  monitor  the  changes  of  different  components  of  the  system  energy.  The  program 
is  implemented  in  two  modifications  which  use  direct  variations  or  those  in  the  Fourier  space  of 
the  initial  functions  components. 

These  program  versions  are  realized  in  the  Object  Pascal  language  in  Delphi  5  Enterprise 
integrated  development  environment.  The  elementary  surfaces  are  calculated  with  the  linearized 
expressions  (see  Deliverable  2).  The  elementary  volumes  are  calculated  by  numerical  integration 
of  the  corresponding  expressions  with  the  help  of  the  standard  method  of  trapezoids. 

In  the  first  version  we  use  the  method  of  local  variations  under  the  integral  condition  of 
volume  constancy  whereas  in  the  second  -  the  method  of  shape  variations  under  the  same 
condition.  Both  versions  use  the  standard  method  of  unconditional  successive  multidimensional 
optimization  in  the  presence  of  local  steric  restrictions.  The  exemplary  screenshots  are  given 
below  (see  Figs.2-4). 
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Fig.l.  Modelling  of  a  planar  bilayer  (note  a  difference  in  xend  in  the  upper  and  lower 
screenshots). 
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Fig.2.  The  beginning  of  the  optimization  procedure  (method  of  local  variations). 
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Fig.3.  The  final  result  (equilibrium  planar  bilayer). 
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Fig.4.  The  course  of  optimization  in  the  Fourier  space. 


The  advantage  of  the  method  of  shape  variation  consists  in  the  guaranteed  absence  of 
punctual  disruptions  of  derivatives  which  destabilize  the  computational  scheme.  Yet,  each 
iteration  requires  full  recalculation  of  the  system  geometry  and  energy  parameters  of  the  whole 
membrane,  essentially  impeding  the  program  execution. 

In  both  methods  the  principle  of  steric  overlap  exclusion  for  lipid  headgroups  situated  in 
two  neighbouring  discretization  points  holds.  To  test  the  presence  of  such  overlap,  the  following 
relationship  for  each  pair  of  neighbouring  discretization  points  i  and  i  + 1  is  verified: 


|  dt>lt 
l^j+i  >  h+ 1 


(14) 


where  d  is  the  distance  from  the  membrane  midline  to  the  point  of  crossing  of  directors  of 
corresponding  lipids,  l  =  b  +  lh  is  the  whole  length  of  the  lipid  (see  Fig. 5).  If  this  relationship 

does  not  hold,  then  the  corresponding  variation  causing  the  steric  overlap  is  rejected. 

The  volume  constancy  condition  is  introduced  with  the  help  of  the  standard  Lagrange 
method  of  undetermined  multipliers.  The  following  term  is  inserted  into  the  free  energy 
functional: 
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fv=A\{v(R)-v0(R)]dR 


where  A  is  the  Lagrange  multiplier  and  v0(R)dR  is  the  infinitesimal  volume  of  the  planar 
unperturbed  membrane.  The  optimization  problem  is  solved  for  a  series  of  values  A  unless  the 
value  of  noKa  3naHCimc  fv  becomes  sufficiently  close  to  zero.  The  search  of  the  corresponding 

value  of  A  is  basically  a  standard  task  of  finding  the  zero  of  a  function  of  one  variable  and  is 
performed  by  bipartitioning. 


Fig.5.  To  condition  (14). 


4.  Conclusions 

The  defined  problem  of  optimization  of  the  membrane  shape  in  the  presence  of  the  CPP 
has  turned  out  to  be,  methodically  and  technically,  far  more  complex  that  it  could  be  expected 
from  the  existing  literature.  Although  the  chossen  approach,  based  on  the  extended  OFM,  in 
principle  allows  one  to  describe  complex  membrane  structures  of  any  shape,  yet  it  has  been 
previously  exploited  for  calculations  of  mainly  small  deviations  from  the  planar  shape  of  the 
bilayer.  In  those  cases  the  corresponding  Euler-Lagrange  equations  can  be  linearized  and  the 
equilibrium  shape  of  the  membrane  can  be  found  comparatively  easily.  The  general  task  of  large 
nonlinear  defonnations  is  intrinsically  complicated  and  specific  for  the  following  reasons: 

(i)  the  Euler-Lagrange  equations  are  very  complex  in  the  general  case; 

(ii)  the  solution  corresponding  to  the  planar  unperturbed  membrane  represents  a  singularity 
of  these  equations.  Consequently,  this  requires  specific  regularization  procedures; 
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(iii)  some  terms  of  these  equations  contain  integrals  with  their  analytical  expressions  having 
critical  points  in  the  parameters  range  of  interest.  This  forces  to  employ  numerical 
computing  of  the  integrals  and  leads  to  far  lower  computational  intensity; 

(iv)  finally,  the  functions  under  optimization  are  fast-varying  and  therefore  strongly  toughen 
the  requirements  to  numerical  algorithms. 

An  alternative  method  of  optimization  -  method  of  local  variations  -  has  also  turned  out 
quite  complex  and  specific  in  realization.  The  corresponding  complexities  are  connected  with 
bifurcational  dependences  of  the  optimal  solution  on  parameters  values  and  also  with  the 
presence  of  multiple  local  minima,  nonlocal  character  of  the  volume  constancy  condition  as  well 
as  with  the  instability  of  the  standard  computational  scheme  with  respect  to  punctual  "defects". 
In  view  of  this,  we  have  decided  to  develop  both  approaches  in  parallel  since  none  of  them  is 
noticeably  preferable. 

The  results  obtained  are  of  preliminary  character  and  indicate  the  necessity  of  further 
development  of  the  optimization  procedure  for  the  free  energy  functional  in  order  to  construct  a 
stable  scheme,  yielding  first  of  all  the  unperturbed/planar  membrane  shape.  Then  the  numerical 
experiment  on  the  (bilayer  membrane  +  CPP)  system  will  go  comparatively  smoothly.  The 
development  of  sseveral  computational  approaches  seems  quite  reasonable  at  the  present  stage. 

The  need  of  some  further  modernization  of  the  extended  OFM  model  in  order  to  obtain 
stable  solutions  for  the  unperturbed  membrane  looks  also  very  likely.  This  will  be  cleared  up  in 
the  course  of  further  application  of  the  above-described  versions  of  the  optimization  schemes. 
Probably,  more  powerful  computational  means  than  presently  used  standard  PCs  will  be  needed. 

Appendix 

Here  we  decipher  the  designations  used  above,  reflecting  the  complexity  of  the  problem. 
Recall  that  the  superscripts  of  monolayers  (supposing  identical  expressions  for  both  under  the 
replacement  b  — >  q  )  and  the  bars  denoting  renormalized  dimensionless  quantities  are  omitted. 
Besides,  we  also  omit  the  terms  of  the  Lagrangian  that  contain  R  only  (as  full  derivatives  with 
respect  to  R  ,  inessential  for  minimization)  and  denote 

s(R;b)  =  s,  s(R;b  +  lh)  =  s. 

The  first  derivatives  of  L  read 
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Besides, 


t\  =  Rz'I2  +  (Rb'z  +  bRz ')  /3  +  bRb'zI4 
t2  =  RI3  +  ( RbR  +  bR )  /4  +  bRbRI5 
h  =  Rz'h  +{Rb'z+  bRz ')  /4  +  bRb%, 


and  of  course  b  =  \Jbl+b:  . 

Below  we  list  the  derivatives  entering  Eqs.(7). 
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The  volume  constancy  condition  requires  derivatives  —  and  also  various  second  partial 

d// 

derivatives  of  v.  Here  they  are: 


dv  ,  dv  ,  .  dv  j,T  ,  T  dv  .  T 
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d2v  _  4  d2v  _  4  d2v  _  4 
dbzdz'  ~  1  db:db'R  ~  2  db:db':  ~  3 


S2v  _h  dh 

dz'2  2  dz' 

=  K \RI2  +  bRI3-  Rz'2J2  -z'(2Rb'z  +  bRz')J3-b'z(  Rb':  +  2bRz')J4-  bRb'2J5  \ 

d2 v  _  dt2  _  b  8tx 
dz'db'R  ~bz!k!~bzWR 

=  bz  [ -Rz'J3  -  ( Rb':  +  Rb'Rz'  +  bRz')J4  -  ( Rb'Rb'z  +  bRb'z  +  bRb'Rz')J5  -  bRb'Rb':J6] 

d2 v  dt3  _  ^  8tx 

dz'dbl  2  dz'  2  db'z 

=  bz  [RI3  +  bRI4  -  Rz,2J 3  -  z'{2Rb'z  +  bRz')J4  -b'z  ( Rb'z  +  2 bRz')J5  -  bRb'2J~_ 


d\  _  d\ 
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see  above 


bz  RI4+bRI5  RJ4  ( 2  RbR  +  bR )  J5  bR  ( RbR  +  2b R  )  J(]  bRbR  J 2 


d2v  dt3  .  dt2 

eb’sd  b[  '*4  -S % 
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see  above 


see  above 


d2v  _  b 
8b'2  2  8b': 

=  bz[R14  +  Vs  -Rz'2j4  ~z'{2Rb'z  +  bRz')J5  -b'z(Rb'z  +  2bRz')J6  - bRb'2J1 
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Auxiliary  relationships: 


~z'Jn  ~  b'zJn+\ 


--J  -til  n  -  -z'T  -til 

^  n+ 1  n+ 2  r  ^  ^  n+ 1  n+ 2 
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etc.  So  that,  e.g., 
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dbR  dbz 
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so  forth. 
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1.  Introduction 

In  this  last  Deliverable  we  present  some  results  of  calculations  within  the  framework  of 
our  modification  of  the  opposing  forces  model  (OFM),  in  particular,  with  allowance  of  tilt 
rigidity  of  lipid  layers.  The  changes  of  the  membrane  shape  in  the  course  of  peptide  binding  and 
further  insertion  into  the  membrane  are  described.  The  calculation  program  allows  one  to 
investigate  the  influence  of  ionic  strength,  peptide  charge  and  shape,  membrane  parameters,  etc, 
upon  the  binding  energy  and  the  system  energy  profile.  At  the  same  time,  close  attention  is 
focused  to  the  revealed  limitations  and  inconsistencies  of  some  points  of  the  OFM  that  impede 
its  full-value  application  to  the  analysis  of  the  translocation  process.  Thus,  apart  from  necessity 
of  involving  the  tilt  rigidity,  the  topological  and  energetic  drawbacks  of  the  OFM  are  shown,  as 
well  as  its  criticality  to  instabilities  of  lipid  aggregates  described  within  its  framework  (first  of 
all,  of  a  planar  (bi)layer)  and  to  non-locality  of  some  of  the  OFM  basic  relationships.  The  ways 
of  eliminating  these  drawbacks  are  pointed  out.  Proceeding  from  all  this,  several  promising  leads 
of  further  research  are  analyzed  in  detail. 

2.  Computational  scheme  and  results 

It  has  become  evident  in  the  course  of  modeling  that  the  extended  OFM  has  several 
serious  drawbacks.  In  particular,  with  the  typically  exploited  free  energy  functional  [1]  it  is 
unable  to  produce  the  unperturbed  shape  of  a  planar  membrane  (see  below).  Apart  from  solving 
these  problems  theoretically,  we  managed  to  achieve  the  convergence  of  computational  scheme 
by  introducing  the  tilt  rigidity  term  of  ~10  kT  into  the  optimized  functional.  The  volume  of  the 
hydrophobic  part  of  the  membrane  was  maintained  fixed  by  applying  a  standard  scheme  with 
Lagrange  multiplier  X  which  value  was  guessed  by  tries- and-errors.  Luckily,  the  value  of  X 
appears  to  be  independent  of  the  electrostatic  shielding  constant  and  membrane  dimensions,  that 
made  the  calculations  feasible  without  re-evaluation  of  X ,  what  is  extremely  costly 
computationally. 

It  is  necessary  to  note  that  the  shape  of  the  membrane  in  our  computations  is  presented  as 
a  series  of  few  cosines;  therefore,  possible  deformations  of  the  membrane  are  limited  to  those 
described  as  a  finite  sum  of  cosines.  Particularly,  point  defects  of  the  lipid  packing  that  require  a 
very  large  number  of  Fourier  components  cannot  emerge  in  calculations.  This  is  a  possible 
reason  of  the  fact  that  the  computational  scheme  is  stable  while  the  underlying  energy  functional 
has  no  stable  steady-states  in  terms  of  Lyapunov  stability  analysis. 
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2.1.  Computational  details 


The  program  that  implements  the  algorithm  of  shape  variations  (see  Deliverable  2,3)  was 
used  with  several  modifications.  The  program  searches  for  five  unknown  functions,  z,  b,  q,  0 
and  c|)  in  the  interval  from  0  to  L.  Each  function  is  expanded  in  a  series  of  M  sines  or  cosines  as 
follows: 


M 

z{R)  =  cos 

1=1 

M 

0(JR)  =  £0~sin(O3,/?) 

1=1 

71  7t(/-l) 

W:  =  -  +— - - ~ 

2  L  L 


where  M  is  the  number  of  harmonics. 

The  same  expansion  as  for  z  is  used  for  b  and  q  and  the  same  expansion  as  for  0  is  used 
for  c|) .  This  choice  of  expansion  coefficients  allows  us  to  satisfy  the  desired  boundary  conditions 


30 
3  R 


R=O.L 


=  0 


z(L)  =  0 


automatically.  The  same  conditions  are  also  applied  to  other  functions.  This  ensures  that  that  the 
studied  segment  of  the  membrane  is  smoothly  connected  to  the  unperturbed  planar  bilayer  at 
point  L  and  there  is  no  tilt  of  lipid  directors  at  point  R= 0  (prohibited  by  symmetry).  The 
derivatives  of  all  functions  are  computed  analytically.  The  program  performs  the  minimization  in 

the  space  of  expansion  coefficients  z,,  bj ,  qt ,  0~  and  (jr  using  the  method  of  simple 

unconditional  minimization  (see  Deliverables  2,3). 

The  electrostatic  interaction  between  the  charged  peptide  and  the  membrane  is  taken  in 
the  form 


N„  L 


''elect 


=11 


1 


QeSu(R) 


rti  4jie0  d  (R)ah(R) 


exp 


dj(R) 


dR  , 


2 


where  Np  is  the  number  of  charges  of  the  peptide,  £0  is  the  vacuum  dielectric  constant,  qe  is  the 
elementary  charge,  su(R)dR  is  the  infinitesimal  area  of  the  interface  at  the  level  of  the  polar  head 
groups,  ah  is  the  area  per  lipid  at  the  same  level,  d0  is  the  shielding  constant  for  electrostatic 
interactions  (similar  but  not  identical  to  the  Debay  length  in  continuous  electrolyte,  and  dj  is  the 
distance  from  the  / - 1 h  charge  of  the  peptide  to  a  given  element  of  the  membrane  surface: 


dj  (R)  =  Jxint(R)2  +  (Z.m  (R)-ZjY 


where  Xint  and  Zint  are  the  coordinates  of  heads  of  the  lipids  anchored  at  point  R;  zj  is  the 

coordinate  of  the  / - th  charge  of  the  peptide  (the  latter  is  positioned  at  a-0).  Multiple  charges  are 
positioned  evenly  along  the  peptide  as  described  below. 

The  close  contact  of  the  peptide  with  the  membrane  is  modeled  by  the  “nearly  hard  core” 
potential  assigned  to  the  peptide: 


E\/dw  (r) 


O'-fcore) 

K. 
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r  >  r 

core 


r  <  r 


+  8 
+  8 


where  rcore  is  the  radius  of  the  hard  sphere  that  surrounds  the  lowermost  charge  of  the  peptide,  k 
and  K  are  adjustable  parameters,  r  is  the  distance  between  the  given  discrete  portion  of  the 
interface  and  the  lowermost  charge  of  the  peptide,  8  is  a  very  small  shift.  In  typical  simulations 
k=  1,  A=1010,  8  =  1 0  8 .  The  common  problem  of  modeling  the  constraints  in  variation  methods  is 
the  “dead  lock”  which  often  occurs  after  the  first  contact  with  the  hard  core.  Our  form  of 
interaction  produces  a  behavior  which  is  almost  indistinguishable  from  that  under  the  true  hard 
core  impermeability  criterion  and  allows  us  to  overcome  computational  complications. 

2.2.  Shape  of  the  membrane 

The  following  series  of  snapshots  illustrates  the  principal  behavior  of  the  membrane  with 

o 

the  peptide  approaching.  The  value  of  the  electrostatic  shielding  constant  is  kept  very  high  (25  A) 
in  order  to  make  the  shape  variations  more  pronounced.  At  the  first  stage  the  membrane  remains 
unaffected  in  its  equilibrium  state  (Fig.  la).  When  the  peptide  goes  closer  to  the  charged  interface 
the  electrostatic  attraction  begins  to  curve  the  membrane  toward  the  peptide  forming  a  bell-like 

shape  maintained  by  the  balance  of  electrostatic  and  elastic  forces  (Fig.  lb).  Eventually  the 
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membrane  collides  with  the  hard  core  of  the  peptide  (Fig.  lc).  Starting  from  this  stage  the  forces 
of  hard-core  van-der-Waals  (VDW)  repulsion  and  electrostatic  attraction  act  in  opposite 
directions.  As  a  result  the  peptide  “digs”  a  central  well  in  the  upward  bell-like  shape  of  the 
membrane  forming  the  "M-iike"  shape  (Fig.  Id).  Since  the  peptide  goes  close  to  the  membrane 
baseline,  the  average  curvature  of  the  M- shaped  membrane  decreases  and  the  associated  elastic 
energy  becomes  small  in  comparison  to  early  bell-like  structures.  However,  the  electrostatic 
interaction  energy  between  the  peptide  and  the  charges  membrane  interface  remains  high.  As  a 
result,  the  M-like  structures  correspond  to  the  minimal  energy  of  the  membrane-peptide  system 
and  are  likely  to  be  observed  in  reality.  Further  approach  of  the  peptide  forces  the  membrane  to 
follow  it.  The  wings  of  the  M-shape  become  flattened  and  the  central  depression  under  the 
peptide  becomes  the  dominant  feature  (Fig.  le).  Further  movements  of  the  peptide  lead  to  the 
deepening  of  this  depression  and  to  dramatic  increase  of  the  elastic  deformation  energy  (Fig.  If). 
The  plots  under  the  schemes  of  the  membrane  show  the  deviations  of  the  lipid  director  length 
(blue  -  the  upper  monolayer,  grey  -  the  lower)  and  the  director  tilt  angles  (red  -  the  upper 
monolayer,  green  -  the  lower).  The  deviations  are  not  scaled  and  serve  for  comparative  purposes 
only. 

o 

The  value  of  the  shielding  constant  d=  25  A  is  not  realistic  since  it  corresponds  to  the 

o 

environment  of  very  low  ionic  strength.  More  realistic  values  are  d=  3-=-6  A.  For  such  values  the 
overall  picture  of  the  membrane  deformations  remains  the  same  but  becomes  much  less 
pronounced. 


0 

a)  No  interaction 


b)  Formation  of  the 
bell-like  shape 
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c)  First  hard-core 
contact 


d)  The  M- shape 


e)  Flattening  of  f)  Depression 

the  M- shape 


Fig.  1.  Snapshots  of  the  membrane  shape  for  various  positions  of  the  peptide  (see  the 

o  o 

text).  Radius  of  the  modeled  part  of  the  membrane  is  140  A  (shown  cross-section  is  280  A  wide). 

o 

The  peptide  is  modeled  by  a  hard  10  A  sphere  with  the  single  point  charge  in  its  center.  Vertical 
scale  is  double- stretched  to  make  the  shape  variations  more  pronounced. 
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2.3.  The  influence  of  ionic  strength 


We  calculated  the  potential  energy  profiles  of  binding  for  different  values  of  the  shielding 
constant.  The  results  are  shown  in  Fig. 2.  It  is  clearly  seen  that  the  shielding  constant  influences 
the  shape  and  the  depth  of  the  energy  profile  dramatically.  Decrease  of  tine  shielding  constant 
decreases  the  binding  energy  in  non-linear  fashion  (Fig. 2).  The  most  physiologically  relevant 
value  of  the  shielding  constant  is  rather  hard  to  estimate  since  it  includes  not  only  the 
conventional  Debay  factor,  but  also  the  effective  influence  of  the  polarization  effects  and  local 

o 

conformational  changes  in  the  membrane  and  the  peptide.  The  value  of  5Acan  be  assigned  as  a 
reasonable  estimate  that  leads  to  the  binding  energy  of  7.7  kT.  This  value  is  quite  close  to 
experimentally  determined  binding  energy  per  peptide.  It  is  necessary  to  emphasize,  however, 
that  this  value  is  obtained  for  a  single  point  charge  approaching  to  the  membrane,  thus  it  can  be 
considered  as  a  rough  estimate  only. 

2.4.  The  influence  of  the  peptide  charge 


To  study  the  influence  of  distributed  peptide  charge  on  the  binding  energy,  we  adopted 

o 

the  following  model  of  the  peptide.  The  latter  is  represented  by  a  vertical  line  of  length  L=15  A. 
The  first  charge  is  located  at  the  lower  end  of  the  line,  while  others  are  distributed  at  even 
distances  along  the  line.  The  hard-core  repulsion  was  calculated  for  the  lowest  charge  only  since 
other  charges  do  not  contact  with  the  membrane  surface  directly.  The  shielding  constant  is  fixed 

o 

at  the  value  of  5 A  It  can  be  seen  in  Fig.  3  that  the  increase  in  the  peptide  charge  leads  to  gradual 
(almost  linear)  increase  in  the  binding  energy. 
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Fig.  2.  Energy  profiles  and  binding  energy  at  different  levels  of  shielding. 
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Fig.  3.  Energy  profiles  and  the  binding  energy  for  various  peptide  charges. 

3.  The  OFM  model  should  be  further  modified  to  describe  pore  formation 

The  ultimate  goal  of  the  current  study  is  to  simulate  the  formation  of  the  pore  in  the 
membrane  caused  by  the  influence  of  the  CPP.  The  most  probable  mechanism  of  pore  formation 
is  formation  of  the  semi-toroidal  micelle-like  rim  around  the  water-filled  pore.  This  excludes 
unfavorable  contacts  of  lipid  tails  with  water  and  forms  a  hydrophilic  channel  for  the  permeating 
peptide.  Formation  of  such  a  membrane  shape  implies  that  the  packing  of  the  lipid  tails  is  not 
uniform  and  differ  considerably  in  the  micellar  rim  and  in  the  distant  parts  of  planar-  bilayer.  It 
was  claimed  in  several  publications  (see  Refs,  in  Deliverables  2,3)  that  the  modified  semi- 
phenomenological  OFM  model  is  applicable  to  membrane  shapes  with  non-uniform  chain 
packing.  Particularly,  the  energy  of  the  semi -micellar  bilayer  edge  was  computed  [1],  However, 
there  is  one  extremely  important  assumption  in  these  works,  not  emphasized  by  the  authors  but 
critical  for  applications  of  the  modified  OFM  model.  In  all  mentioned  works  the  existence  of  a 
stable  symmetric  bilayer  is  postulated,  not  derived  from  the  first  principles.  This  approach  looks 
sufficient  if  the  geometry  of  the  membrane  system  under  study  is  simple  enough  and  known  from 
experiment.  However,  it  fails  if  there  is  no  a  priori  knowledge  of  the  membrane  shape.  The 
binding  problem  implies  that  the  membrane  remains  in  the  bilayer  form  but  bends  to 
accommodate  the  peptide.  This  process  can  be  classified  as  manageable  by  the  OFM- like  model 
we  use.  However,  the  problem  of  pore  formation  is  an  unfavorable  case  when  there  is  no  notion 
of  the  shape  of  intermediate  states  of  the  membrane. 

Fet  us  examine  the  problem  in  more  details.  The  transition  from  a  planar  bilayer  to  a 
semi-toroidal  pore  changes  the  topology  of  the  membrane  surface  qualitatively.  In  rigirous 
topological  terms ,  a  planar  membrane  is  isomorphous  to  a  sphere,  while  a  pore  is  isomorphous  to 
a  toroid.  Within  the  OFM  the  membrane  is  treated  in  continuous  approximation  and  the  lipid 
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directors  are  “anchored”  on  the  membrane  midline.  However,  the  latter  becomes  discontinuous 
when  the  membrane  evolves  from  a  bilayer  to  a  pore.  Since  all  energy  contributions  are 
expressed  in  terms  of  the  midline,  they  are  not  topologically  invariant.  The  simplest  illustration  is 
given  in  Fig.  4a.  In  the  “ideal”  membrane  edge  capped  by  the  micellar  region  the  directors  of  all 
lipids  that  form  the  micellar  cap  originate  from  the  same  point  (the  origin  of  the  membrane 
midline).  At  the  same  time,  in  the  bilayer  part,  only  one  lipid  per  monolayer  can  originate  from 
the  same  point.  In  Ref.  1  this  problem  was  solved  by  assuming  that  the  micellar  cap  is  already 
formed.  Then  the  system  was  split  into  two  topologically  different  parts  -the  bilayer  and  the 
micelle  and  different  energy  functionals  were  written  for  each  of  them.  However,  we  can  not  go 
this  way  because  we  have  to  describe  the  formation  of  the  micellar  part,  that  is,  a  smooth 
transition  from  the  geometry  of  a  bilayer  to  the  geometry  of  a  micelle.  The  only  possible  way  to 
do  this  is  illustrated  in  Fig.  4.  If  the  number  of  discrete  points  that  represent  the  membrane  is 
large  enough,  then  the  difference  between  these  two  representations  becomes  negligible. 

However,  it  is  still  impossible  to  describe  a  smooth  transition  from  the  bilayer  to  the  pore. 
In  order  to  form  the  “seed”  of  the  pore,  the  origin  of  the  midline  (bold  dot  in  Fig.4)  should  move 
away  from  zero  (symmetry  axis  of  the  future  pore).  In  this  case  there  is  no  lipid  “anchored” 
between  the  pore  axis  and  the  midline  origin.  From  the  physical  point  of  view  this  means  that  the 
tails  of  the  lipids  closest  to  the  axis  will  repack  themselves  to  fill  the  gap.  Since  this  packing  will 
obviously  be  not  ideal,  the  system  energy  will  increase.  However,  appearance  of  such  a  gap 
becomes  fatal  for  the  OFM.  The  energy  of  the  membrane  is  computed  as  an  integral  over  the 
midline,  and  the  regions  where  the  latter  is  absent  contribute  nothing  to  the  total  energy.  In  other 
words,  the  “missed  volume”  shown  in  Fig.4  has  zero  energy  in  terms  of  the  OFM.  The  system  is 
absolutely  invariant  to  the  size  and  shape  of  this  region,  what  contradicts  physical  reality. 

The  next  step  of  pore  formation  is  the  rapture  of  the  interfaces.  This  allows  water  to  fill 
the  pore  and  to  contact  with  the  hydrophobic  tails  of  the  lipids.  In  reality,  such  a  contact  is 
extremely  unfavorable  and  the  first  lipids  of  the  upper  and  lower  interfaces  will  approach  each 
other  to  eliminate  the  hydrophobic  mismatch.  However,  the  OFM  contains  no  explicit  terms 
describing  the  hydrophobic  mismatch  -  the  system  does  not  “feel”  that  the  tails  are  exposed  to 
water. 

There  is  one  more  energy  term  missed  in  the  OFM.  The  heads  of  the  lipids  are  often 
charged  and  repel  each  other.  This  interaction  within  single  interface  is  implicitly  described  by 
the  B  /  ah  term  of  the  OFM.  However,  if  two  separate  interfaces  are  located  close  to  each  other, 

there  is  a  strong  electrostatic  interaction  between  them.  This  interaction  should  be  described 
explicitly  but  it  is  missed  in  the  OFM.  Introduction  of  this  additional  energy  term  is  far  from 
trivial.  The  micellar  cap  formally  consists  of  two  parts  -  one  from  the  upper  and  another  from  the 
lower  interface.  An  explicit  electrostatic  interaction  should  act  between  these  parts,  but  at  the 
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same  time  the  interface  of  the  micelle  is  topologically  continuous  and  should  be  described  by  the 
implicit  term  itself. 

All  the  three  energy  components  missed  in  the  OFM  are  schematically  shown  in  Fig.4. 

4.  Internal  inconsistencies  of  the  extended  OFM 

4.1.  The  tilt  energy 

It  is  well  known  that  lipid  bilayers  possess  some  rigidity  to  tilt  deformations  that  are  often 
described  in  the  linear  approximation  by  the  tilt  modulus  (see  e.g.  [2]).  The  physical  meaning  of 
tilt  rigidity  in  terms  of  the  OFM  originates  from  the  fact  that  the  fluctuating  hydrophobic  tails  of 
the  lipids  cannot  cross  the  interface  with  water  in  order  not  to  be  exposed.  If  the  director  of  a 
lipid  is  tilted  from  the  normal  of  the  hydrophobic-water  interface,  a  certain  part  of  the  space 
available  for  the  fluctuating  hydrocarbon  chain  in  normal  orientation  becomes  excluded.  This 
leads  to  some  loss  of  conformational  freedom  and,  as  a  result,  to  the  change  in  the  entropic 
contribution  to  the  free  energy  of  the  lipid.  The  energy  term  associated  with  this  effect  can  be 
written  approximately  as  Emt  =  ktilia 2 ,  where  a  is  the  angle  between  the  normal  to  the 

hydrophobic-water  interface  and  the  lipid  director,  and  ktn,  is  the  tilt  modulus.  It  was  shown  that 
the  tilt  rigidity  of  the  bilayers  is  quite  high  and  the  tilt  modulus  reaches  -10  kT  [2].  However,  this 
term  was  not  included  into  the  modified  OFM.  We  have  found  that  this  term  is  vital  for 
maintaining  the  stability  of  the  membrane.  If  this  term  is  omitted,  certain  “unfortunate”  situations 
can  appear  during  the  optimization  when  the  lipids  orient  themselves  almost  parallel  to  the 
interface  and  get  stuck  in  this  geometry  (see  Fig. 5).  This  configuration  is  clearly  nonphysical, 
since  the  hydrocarbon  tails  have  significant  thickness  and  will  be  exposed  to  water  if  a  is  close 
to  90°.  We  introduced  the  tilt  term  into  our  simulations  and  assigned  the  value  10  kT  to  the  tilt 
modulus. 

4.2.  The  challenge  of  non-locality 

In  order  to  find  the  area  and  the  volume  per  lipid  the  following  expressions  are  used  (for  the  sake 
of  clarity  we  here  consider  a  quasi-one-dimensional  case  that  does  not  restrict  the  problem 
generality): 


9 


Unperturbed  bilayer 


Origin  of  the  midline 
moves  from  zero,  forming 
a  seed  of  micellar 
topology,  but  the  interfaces 
remain  continuous. 


Origin  of  the  midline 
moves  further,  the 
interfaces  break,  the 
semitoroidal  rim  forms,  but 
the  tails  of  lipids  at  the  edge 
are  exposed  to  water. 


Former  upper  and  lower 
interfaces  merge  forming  a 
complete  semi -toroidal 
micelle-like  rim  with  no 
hydrophobic  mismatch. 


Fig.  4. 
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Excluded  volume 


Fig.  5. 


R+AR 

at(R)  =  |  s(Rr)dR' »  s(R)AR 

R 

R+AR 

D=  J  v(Rr)dR'^v(R)AR 

R 

what  is  approximately  valid  if  A R  «  Rmin ,  where  Rmin  is  the  minimal  local  curvature  radius  of 
the  monolayer  in  the  course  of  membrane  deformation.  For  large  deformations  this 
approximation  is  not  valid  and,  strictly  speaking,  both  equations  should  be  solved  to  find  A R  . 
The  latter  is  the  radial  size  of  an  average  lipid  anchored  at  point  R  of  the  midline.  It  is  easy  to 
show  that  for  an  unperturbed  bilayer  A R  ~  b0 ,  where  bo  is  the  equilibrium  thickness  of  a 

o  o 

monolayer  (typically  ~  10  A).  This  means  that  an  object  that  spans  the  radial  dimension  of  ~10  A 
is  treated  as  an  infinitesimally  small  volume  which  is  integrated  over  R  to  obtain  the  total  energy 
of  the  monolayer  FM  (see  e.g.  Eq.(12)  in  Deliverable  2). 

However,  it  has  turned  out  to  be  a  rather  rough  approximation.  In  fact,  in  the  micellar  part 
of  the  system  A R  ~  Rmw  and  the  area  per  lipid  should  be  expressed  by  the  functional 

R+AR(R) 

a^R)  =  J  s (R')dR' .  In  this  case  the  minimization  of  the  functional  FM  can  no  longer  be 

R 

performed  by  the  standard  Euler-Lagrange  method.  The  reason  is  that  the  later  can  only  be 
applied  if  the  functional  of  form  F(x,t )  =  /'(x,x  t)clt  is  minimized  while  in  our  case  we  have  a 
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“functional  of  a  functional”.  Minimization  of  such  a  complex  construct  is  itself  a  non-standard 
and  challenging  mathematical  problem.  Additional  complication  comes  from  the  fact  the 
computational  methods  of  optimization  of  such  complex  “double”  functional  are  underdeveloped 
(our  results  to  these  specific  issues  are  now  in  progress  and  will  be  presented  elsewhere). 

4.3.  The  ambiguity  of  the  constant  volume  constraint 

A  physical  formulation  of  our  problem  implies  the  “flat  unperturbed  membrane”  as  a 
reference  state.  This  trivial  phrase  contains  a  subtle  pitfall,  however.  A  real  biological  membrane 
(say,  a  liposome  of  large  radius)  has  a  fixed  total  volume  of  the  hydrophobic  core,  conserved 
with  high  accuracy.  We  consider  a  much  smaller  piece  of  the  membrane  and  imply  that  it 
remains  unperturbed  infinitely  far  from  the  point  of  perturbation  by  the  peptide.  In  practice,  we 

o 

consider  the  fraction  of  the  membrane  spanning  from  -/  to  /,  where  l~  100- 200  A  In  the 
unperturbed  state  this  piece  contains  N  lipids  and  has  volume  V  of  the  hydrophobic  core. 
Perturbation  of  the  membrane  causes  the  changes  in  chain  packing  and  redistribution  of  lipids 
along  the  membrane.  Suppose  that  the  membrane  takes  a  bell-like  shape  (see  Fig. 6).  The  lipids 
located  at  points  /  and  -/  tend  to  move  toward  the  center  of  perturbation.  Several  kinds  of  the 
boundary  conditions  are  then  possible: 

(i)  Impermeable  walls  -  no  lipids  can  cross  the  boundaries,  /  =  const ,  N  =  const , 

V  =  const .  If  the  membrane  is  deformed,  its  area  increases  and  substantial  surface  tension 
appears.  The  energy  of  this  artificial  strain  is  added  to  the  “intrinsic”  deformation  energy  of  the 
membrane  introducing  some  systematic  error. 

(ii)  “Transparent”  walls  -  the  lipids  can  freely  cross  the  boundaries,  /  =  const  but  N  and 

V  are  not  conserved.  This  seems  to  be  the  most  physically  correct  version  -  the  lipids  can  leave 
the  given  volume  but  accumulate  outside  in  some  other  parts  of  a  huge  liposome.  In  practice,  this 
leads  to  the  fact  that  all  lipids  tend  to  leave  the  volume.  The  reason  of  this  is  that  the  energies  in 
the  OFM  are  always  positive  and,  therefore,  the  state  with  no  lipid  is,  erroneously,  the  global 
energy  minimum.  Thus,  this  type  of  boundary  conditions  is  not  applicable. 

(iii)  Moving  impermeable  walls  -  no  lipids  can  cross  the  boundaries,  N  =  const , 

V  =  const ,  but  the  value  of  /  can  be  adjusted  to  reduce  the  tension  caused  by  deformation.  This 
type  of  boundary  conditions  is  the  most  preferable  since  it  retains  the  integrity  of  the  membrane 
and  eliminates  artificial  strain.  Unfortunately,  moving  boundaries  are  quite  hard  to  implement 
computationally  because  the  discrete  lattice  used  in  simulations  should  be  redefined  at  each 
iteration.  This  results  in  numerous  specific  problems  and  greatly  increases  computational  burden. 
That  is  why  we  decided  not  to  implement  this  method  until  we  improve  the  accuracy  of  the 
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membrane  model  to  the  point  when  the  error  introduced  by  additional  mechanical  strain  will 
become  detectable. 


Impermeable 

walls 

Additional 
mechanical 
strain  appears 


Transparent 

walls 

Planar  membrane  is 
unstable  -  the  number 
of  lipids  decreases 
because  the  energy 
minimum  corresponds 
to  the  absence  of  lipids 


Fig.  6. 


4.4.  Stability  problem 


In  this  subsection  we  elucidate  a  crucial  point  at  issue  for  the  OFM  validity,  namely,  the 
stability  of  lipid  aggregates  within  the  OFM  framework. 

To  illustrate  the  problems  arising  in  the  procedure  of  minimization  of  the  free  energy 

functional,  we  first  consider  a  Lagrangian  L  (y(x),y'(x))  entering  functional  j  L  (y,  y)dx  and 
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provided  that  L  does  not  depend  on  x  explicitly.  Obviously,  the  "stationary"  solution  ys  of 
Eq.(l)  is  defined  by  the  equation 


dl_ 

dy 


=  0. 
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Here  and  below  the  subscript  s  labels  the  functions  taken  at  y  =yJ,y'  =  0,  that  is,  for 
(|)(  v,  y')  4>v  =  (|)  C y  v  - 0 ) .  In  the  standard  linear  analysis  of  stability  one  re-writes  (1)  as  a  set  of  two 
first-order  equations  for  variables  yl  =  y,  y2  =  y  '■ 
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u  dF2 

so  that 


3^2 


■0  because 
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3vj 


■0 ,  see  Eq.(2).  This  means  nothing  else  than  stationary  states  of 


the  variational  problem  with  such  kind  of  Lagrangian  are  always  unstable  (as  it  should  be 
expected  from  a  corresponding  mechanical  analogue).  In  our  context,  we  thus  conclude  that  the 
OFM  is  not  able  to  ensure  the  existence  of  a  stable  planar  (bi)layer! 

One  of  the  ways  out  of  the  situation  is  to  insert  an  artificial  decay  term  proportional  to 
(-/)  in  the  r.h.s.  of  (1),  exactly  as  the  friction  term  is  introduced  in  equations  of  classical 

mechanics.  Of  course,  in  our  case  physical  meaning  could  hardly  be  attached  to  such  "friction".  It 
is  worth  to  note,  however,  that  this  term  affects  neither  the  stationary  solution  ys  nor  the  stability 
domain  in  the  parameter  space  (in  particular,  the  coefficient  P  >  0  in  this  formal  term  -p  y  can 
be  even  infinitesimal).  Indeed,  if  we  suppose  that  now 
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Return  now  to  the  simplest  case  of  a  lipid  layer  within  the  extended  OFM.  To  analyze  the 
problem,  there  is  no  need  in  all  the  five  sought  functions  of  our  model  of  the  membrane.  We  can 
come  to  relevant  conclusions  by  considering  a  quasi-one-dimensional  model  of  a  lipid  layer  in 
which  all  the  lipids  are  perpendicular  to  the  baseline,  and  the  latter  is  fixed  straight.  Then  the 
only  variable  under  consideration  is  the  lipid  length  b  (so  that  y  =  b,  y'  =  b').  In  the  reduced 
(dimensionless)  units  introduced  before  the  corresponding  functional  reads 
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where  the  second  integral  represents  the  lipid  volume  constancy  condition  bdx  -V  ,  X  is  the 

Lagrange  multiplier,  |i  =  (1  -  B)/(l  -  lc) ,  and  1(b)  is  the  layer  length  ("movable  walls");  here  we 

omit  the  bars  over  dimensionless  parameters,  see  Deliverable  2.  The  Euler-Lagrange  equation 
reads 
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and  its  stationary  solution  is  defined  by  the  equation  X  +  — 

db 


0 ,  i.e.  by 


b'=0 


-j^-2  b,(B-tllc)-^  =  X 


(5) 


which  can  be  also  regarded  as  an  equation  to  determine  X ,  for  in  this  case  bs  can  be  found  from 

rl(b) 

the  following  simple  considerations.  Precisely,  for  the  planar  layer  the  condition  bdx  =  V 

rl(b ) 

immediately  results  in  the  relationship  l(bs)  =  V  /  bs.  Recall  now  that  Ldx  is  nothing  else 

J  0 

than  f  } —bdx  .If  b  =  bs ,  then  this  integral  is  reduced  simply  to  —bsl  (bs )  =  —V  =  Nf0  where 
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N  is  the  number  of  lipids.  Therefore,  for  a  planar  layer  the  minimum  of  the  functional 
corresponds  to  the  minimum  of  /0 ,  that  is  bs  =  b0 ;  in  our  dimensionless  units  b0  =  1  (see 


Deliverable  2).  In  particular,  from  (5)  one  has  that  X  -  -2 B  + 


The  stability  condition  (4)  now  reads 
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As  Lhh(  1,0)  =  — - — - — £  and  L^(1,0)  =  1  -B  ,  we  have 
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3-5-2/ 


(6) 


(l-/c)(l-5) 
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Fig.7  shows  the  quantity  (c)F2  / r)_v ) |  as  a  function  of  parameters  B,lc  and  the 
corresponding  stability  domain  defined  by  (6). 


Fig.  7. 


The  return  of  the  parameters  to  their  initial  dimensions  is  not  trivial  because  bQ  itself  is 

defined  as  a  root  of  the  cubic  equation  Txv/y’  +  (B-2ivlc)b^ -yv2  =0  (see  Deliverable  2, 

Eq.(13));  here  the  coefficients  are  bar-less  again.  Of  course,  this  conversion  can  be  performed 
numerically.  An  example  for  the  parameter  values  typically  used  for  lipid  membranes  is  shown  in 
Fig. 8.  One  can  see  how  delicate  the  problem  of  stability  of  the  planar  layer  is. 


B 


Fig.  8. 

x  =  0. 1  kTA A2;  y  =  0.12  JcT/A2;  v  =  432  A3. 


17 


5.  Research  perspectives 


The  present  pilot  project  allowed  us  to  evaluate  one  of  the  possible  approaches  to 
description  of  the  CPP  permeation.  This  approach  is  based  on  the  modified  semi- 
phenomenological  OFM-like  model.  This  model,  to  our  knowledge,  is  the  most  detailed  among 
available  continuous  membrane  models  and  pretends  to  describe  non-uniform  chain  packing  of 
various  membrane  topologies.  However,  we  revealed  fundamental  obstacles  which  discouraged 
us  from  using  this  model  "as  is",  without  substantial  modification  or  even  development  of  new 
concepts.  We  see  the  following  three  directions  of  research  that  can  be  further  developed 
proceeding  from  the  results  of  the  project: 

A.  Development  of  a  continuous  semi -phenomenological  model,  free  from  the  drawbacks 
of  OFM-like  models. 

Requirements'. 

1.  The  energy  function  contains  both  surface  and  volume  terms  sensitive  to  the  local 
geometry. 

2.  The  lipid  directors  are  anchored  on  the  hydrophobic-polar  surface  that  remains 
continuous  in  all  possible  membrane  perturbations  and  rearrangements. 

3.  The  volume  energy  terms  allow  for  possible  not  ideal  packing  of  lipid  tails  and  for 
interactions  between  different  monolayers. 

4.  The  relationship  between  the  area  and  volume  of  the  average  lipid  anchored  in  the  given 
point  is  non-local.  It  depends  on  the  geometry  of  the  certain  membrane  segment  spanned 
by  the  lipid,  not  only  on  the  geometry  at  the  anchoring  point. 

5.  Additional  surface  terms  describing  long-range  electrostatic  interactions  should  be  added. 

6.  A  flat  unperturbed  monolayer  is  a  stable  steady- state  without  any  additional  constraints  or 
assumptions. 

Benefits: 

1.  Simple  enough  to  find  the  optimal  membrane  shape  of  a  perturbed  membrane  by  low-cost 
computational  techniques. 

2.  The  variety  of  parameters  that  correspond  to  real  membranes  and  CPPs  can  be  testes. 
Limitations: 

1.  It  is  not  clear  if  all  the  requirements  can  be  met  simultaneously. 

2.  The  results  obtained  are  only  semi -quantitative  and  depend  strongly  on  the  quality  of 
empirical  parameters. 

3.  Hard  to  correlate  some  empirical  parameters  with  experimental  data. 

4.  Can  never  answer  the  question  about  atomic-scale  details  of  permeation. 
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B.  Development  of  a  highly  simplified  discrete  molecular- level  model  of  the  membrane. 

Requirements'. 

1.  Each  lipid  is  modeled  as  a  simple  rod-like  object. 

2.  An  empirical  force  field  should  be  developed  to  describe  the  interactions  between  the 
rods  in  order  to  mimic  the  basic  behavior  of  lipids. 

3.  The  model  should  mimic  only  the  very  basic  mechanical  properties  of  the  membranes,  no 
need  to  describe  smaller  details. 

4.  The  charges  of  the  polar  heads  are  included  explicitly  as  well  as  electrostatic  interaction 
with  the  peptide. 

Benefits: 

1.  Eliminates  the  topological  problems  of  the  OFM-like  models  because  the  lipids  are 
modeled  at  the  molecular  level. 

2.  Can  model  any  shape  of  the  membrane,  provided  that  a  correct  force  field  is  supplied. 

3.  Specific  interactions  of  the  peptide  with  the  membrane  (like  strong  binding  of  arginines  to 
phosphates)  can  be  taken  into  account  by  additional  energy  terms. 

4.  Can  be  used  to  model  the  permeation  in  dynamics,  not  just  to  evaluate  the  energy  profiles 
of  permeation. 

5.  It  is  much  cheaper  computationally  than  molecular  dynamics,  allowing  longer  simulations 
and  experiments  with  various  parameters. 

Limitations : 

1.  The  development  and  validation  of  the  force  field  is  a  quite  tedious  task. 

2.  The  force  fields  for  rod-like  particles  are  inconsistent  with  existent  molecular  simulation 
software.  Modification  of  the  existing  programs  or  cle  novo  design  of  the  simulation 
package  is  required. 

3.  The  results  obtained  are  still  semi -quantitative  only  and  depend  strongly  on  the  quality  of 
the  force-field. 

C.  Molecular  simulations  at  the  atomic  or  nearly-atomic  level. 

Requirements'. 

1.  The  methodology  of  such  simulations  is  well  known  in  general  but  should  be  adapted  to  a 
particular  task. 

2.  The  protocol  of  pulling  the  peptide  through  the  membrane  should  be  developed  because 
spontaneous  permeation  is  beyond  the  reachable  timescale. 

3.  The  method  should  be  used  in  conjunction  with  analytical  techniques  to  analyze  the  free 
energy  profiles  of  permeation  and  to  calculate  the  transition  rates. 
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Benefits: 

1.  Well-defined  methodology  and  extensively  tested  force-fields  for  the  membrane  and  the 
peptide. 

2.  Atomic-level  details  of  all  binding  and  permeation  events. 

3.  Explicit  comparative  tests  of  real  peptides  and  membranes  are  possible. 

4.  We  have  established  contacts  with  a  world-leading  molecular  dynamics  group  in 
Groningen.  It  is  possible  to  negotiate  the  usage  of  their  experience  and  technical 
capabilities. 

Limitations: 

1.  The  method  of  molecular  dynamics  is  extremely  demanding  computationally  in 
comparison  with  previous  two  ones.  Supercomputers/clusters  are  needed. 

2.  The  time  scale  of  simulations  is  limited  to  ~100ns  what  is  insufficient  to  observe 
spontaneous  permeation  events.  The  “pulling”  should  be  used,  allowing  to  produce  the 
approximate  free  energy  profiles  of  permeation  only. 
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Summary 


The  main  purpose  of  the  present  project  was  a  theoretical  investigation  of  the  possibility 
of  translocation  of  some  positively  charged  peptides  (so  called  cell- penetrating  peptides,  CPPs) 
through  artificial  phospholipid  bilayer  membranes  and  evaluation  of  the  influence  of  different 
physicochemical  factors  upon  this  process. 

As  a  result  of  analysis  of  the  current,  often  contradictory  experimental  data  it  has  become 
clear  that  reliable  conclusions  on  the  feasibility  of  the  given  non-standard  phenomenon  can  be 
obtained  only  on  the  basis  of  its  physical  modeling  within  a  microscopic  (to  reasonable  extent) 
approach.  Proceeding  from  available  computational  means,  we  have  chosen  the  known  (and 
widely  used  in  current  literature  on  membrane  systems)  opposing  forces  model  which  is 
considered  as  the  most  suitable  for  describing  large  deformations  of  lipid  aggregates  and  their 
complexes  with  proteins. 

We  have  developed  an  OFM  version  aimed  at  a  unified  and  continuous  description  of  the 
process  of  binding  and  translocation  of  the  peptide,  including  complex  changes  of  the  system 
geometry.  In  fact,  a  new  approach  to  the  minimization  of  the  free  energy  functional  of  systems 
with  essential  non-locality  and  a  corresponding  computational  formalism  is  constructed. 
Computer  programs  are  created  which  allows  one  to  elucidate  the  shape  of  a  peptide- membrane 
complex  and  to  calculate  the  binding  energy  and  energy  profile  of  the  system  under  further 
insertion  of  the  peptide  into  the  membrane,  to  investigate  the  dependences  of  these  characteristics 
on  the  ionic  strength,  peptide  charge,  membrane  parameters,  etc. 

At  the  same  time,  in  the  course  of  project  execution  we  have  revealed  a  series  of  essential 
drawbacks  and  inconsistencies  of  the  OFM  which  hinder  its  application  for  a  full  translocation 
process,  especially  at  the  pore  formation  stage.  We  present  a  detailed  analysis  of  these  drawbacks 
and  formulate  the  principles  of  construction  of  a  distinctly  improved  phenomenological  model, 
capable  of  correct  description  of  membrane  rearrangement  topology.  On  this  basis  we  develop 
the  prospective  trends  of  further  research,  rating  them  on  their  success  criteria. 
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